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Introduction 


The  realization  of  structures  that  do  not  scatter  electromagnetic  field,  i.e.  structures  that 
appear  invisible  for  EM  waves  in  outer  space,  is  not  a  new  concept.  The  possibility  of  a  plane 
wave  passing  without  distortions  through  a  structure  with  anisotropic  filling  was  theoretically 
first  investigated  in  1960s,  see  [1].  The  basis  of  the  work  was  the  invariance  property  of 
Maxwell's  equations  with  respect  to  transformation  of  space  metric  and  permeability  and 
permittivity  tensors  of  the  medium.  In  [2]-[5]  it  was  shown  that  for  certain  combinations  of 
permittivities  in  a  two-layer  dielectric  ellipsoid,  the  scattered  field  is  zero.  The  analysis 
revealed  that  at  least  one  of  the  layers  should  have  relative  permittivity  smaller  than  one  (i.e. 
local  negative  polarizability,  which  is  inherent  to  plasmonic  materials  and  plasma-like  ENG 
metamaterials),  and  moreover  such  structures  have  the  advantage  of  not  requiring  anisotropic 
materials.  In  [6],  hard  surfaces  were  used  in  the  design  of  supporting  struts  of  reflector 
antenna  feeds,  to  reduce  their  blockage  in  reflector  systems.  In  this  way  it  is  possible  to 
considerably  reduce  the  scattered  field  for  one  angle  of  incidence,  for  both  polarizations  of 
the  incident  wave.  Several  other  concepts  for  obtaining  invisible  scatterers  were  proposed, 
like  minimum  scattering  antennas  and  active  scatterers  ([7],  [8]). 

The  possibility  of  cloaking  objects  using  a  metamaterial  cover  has  extensively  been 
studied  (see  e.g.  [9]-[13]).  Metamaterials  are  new  composite  materials  that  recently  have 
found  applications  in  the  design  of  antennas  and  microwave  components.  These  special 
materials  can  be  formed  by  periodic  arrangements  of  many  sub-wavelength  inclusions  in  a 
dielectric  environment,  in  such  a  way  as  to  achieve  macroscopic  electromagnetic  or  optical 
properties  that  cannot  be  found  in  nature.  In  the  metamaterial  cloak  approach,  material  has 
been  used  to  render  a  volume  effectively  invisible  to  incident  radiation,  i.e.  to  squeeze  space 
from  a  volume  into  a  shell  surrounding  the  concealment  volume.  The  perfect  cloak  ensures 
that  for  any  incident  field  the  electromagnetic  scattered  field  vanishes  in  the  free-space 
external  to  the  cloaking  shell,  and  the  total  field  vanishes  inside  the  free-space  cavity  of  the 
shell.  Thus,  any  object  placed  in  the  cavity  does  not  perturb  the  electromagnetic  field  outside 
the  cloak,  and  an  external  observer  does  not  see  the  object  and  the  cloak,  like  they  are  absent. 
However,  the  coordinate  transformations  that  are  used  for  cloak  design  do  not  affect  the  form 
of  Maxwell’s  equations,  but  they  influence  permittivity  and  permeability  tensors. 
Consequently,  the  materials  needed  for  building  the  cloak  are  fully  anisotropic  with  spatially- 
varying  constitutive  parameters. 
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The  first  practical  realization  of  a  cloak  was  made  by  Schurig  et  al.  in  2006  [14].  The 
realized  cloak  was  not  fully  anisotropic,  but  it  was  designed  to  work  for  only  one  polarization 
(TMZ  polarization).  By  this  design  of  metamaterial  layers  with  only  the  radial  dependence  of 
the  radial  permeability  component  was  needed,  and  it  was  realized  using  the  split-ring 
resonators.  The  described  cloak  was  prototyped  and  embedded  into  scattering  chamber.  The 
distribution  of  the  vertical  component  of  the  electric  field  around  the  cloak  was  scanned,  and 
the  obtained  results  showed  some  decrease  of  the  scattered  field.  Although  this  experiment 
proved  the  basic  idea  of  cloaking,  the  level  of  scattered  field  has  not  been  quantified. 

Previous  work  in  invisibility,  therefore,  gives  rise  to  several  questions: 

1.  Why  the  developed  cloak  with  the  reduced  variation  of  constitutive  parameters  (i.e. 
cloaks  made  from  uniaxial  materials)  scatters  the  electromagnetic  field,  even  in  the 
ideal  case? 

2.  Is  it  possible  to  realize  a  cloak  with  the  reduced  variation  of  constitutive  parameters 
that  is  entirely  invisible  (at  least  at  the  central  frequency)? 

3.  What  is  the  bandwidth  of  a  cloak  that  is  built  using  metamaterial  layers?  Is  it  possible 
to  enlarge  the  bandwidth  of  such  a  cloak? 

4.  Is  it  possible  to  realize  a  cloak  with  the  reduced  variation  of  constitutive  parameters 
that  work  for  oblique  direction  of  incident  wave? 

5.  Is  it  possible  to  realize  a  cloak  that  is  invisible  to  pulse  excitation  (i.e.  to  radar)? 

6.  Is  it  possible  to  extend  the  single-curved  models  (cloaks)  to  double-curved  structures? 

The  purpose  of  the  project  is,  therefore,  to  provide  further  analysis  of  this  issue. 


The  possible  application  of  the  developed  cloak  will  be  reduction  of  the  scattered  field 
from  some  mechanical  structure.  In  more  details,  there  are  many  situations  where 
electromagnetic  waves  are  obstructed  by  some  mechanical  structure.  The  obstruction  may 
represent  aperture  blockage  causing  increased  sidelobes  and  reduced  gain,  if  the  structure  is  a 
part  of  the  antenna  system  or  if  it  is  placed  close  to  the  antenna.  Examples  of  such  structures 
are  objects  placed  in  the  vicinity  of  the  antenna  system  (e.g.  musts  on  ships  near  radars  or 
communication  systems),  supporting  struts  in  large  reflector  systems,  etc.  In  all  such  systems 
the  frequency,  polarization  and  angle  of  incidence  are  known  which  allows  us  to  construct  a 
special  structure  that  will  reduce  the  scattered  field. 
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2  PROJECT  OBJECTIVE 
AND  REALIZED  OUTCOMES 
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Project  objective  and  realized  outcomes 


Uniaxial  cylindrical  cloaks  have  recently  been  proposed  to  prevent  scattering  of 
electromagnetic  waves,  i.e.  to  render  objects  invisible.  The  proposed  cloaks  with  reduced 
variation  of  constitutive  parameters  suffer  from  nonzero  reflectance,  i.e.  they  are  only  partly 
invisible.  Therefore,  the  purpose  of  a  12-month  effort  is  to  develop  an  analysis  method  and 
the  corresponding  software  for  analyzing  cylinders  made  from  metamaterial-based  concentric 
layers  with  an  application  to  invisible  cloak  realization.  The  main  parameter  of  interest  is  the 
scattered  field,  i.e.  if  the  considered  structure  really  appears  invisible  to  the  incident  wave. 
The  considered  structures  are  modeled  as  a  multilayer  anisotropic  cylinder  where  both 
permittivity  and  permeability  tensors  are  diagonal  in  the  cylindrical  coordinate  system.  The 
active  cloaks  can  be  also  analyzed  in  this  way  since  their  homogeneous  electromagnetic 
model  is  based  on  complex  constitutive  parameters.  We  have  considered  incident 
electromagnetic  waves  coming  both  from  normal  and  oblique  incident  direction.  The 
proposed  analysis  method  is  connected  with  a  suitable  global  optimization  routine,  and  the 
developed  software  package  is  used  to  help  in  cloak  design.  A  possible  application  of  the 
considered  prototype  is  reduction  of  scattered  field  from  mechanical  structures  that  need  to  be 
placed  in  front  of  big  antennas,  such  as  radar  antennas  or  reflector  systems. 

The  work  in  the  project  has  been  divided  in  the  following  tasks: 

1.  Making  a  detailed  investigation  of  properties  of  cloaks  that  are  made  from  uniaxial 
materials,  e.g.  from  metamaterial  layers.  The  investigation  includes  characterization  of 
already  developed  cloaks  and  extensive  investigation  of  the  cloak  bandwidth. 

2.  Extension  of  the  program  for  analyzing  uniaxial  multilayer  cylinders  -  with  the 
extended  version  it  will  be  possible  to  analyze  the  case  when  the  incident  field  has 
oblique  direction  of  propagation  and  oblique  polarization. 

3.  Merging  the  program  for  analyzing  uniaxial  cylindrical  structures  with  the  global 
optimization  program.  As  a  global  optimization  algorithm  we  have  chosen  the  Particle 
Swarm  Optimization  (PSO)  algorithm.  This  is  an  evolutionary  algorithm  similar  to  the 
genetic  algorithm  and  to  the  simulated  annealing,  but  it  operates  on  a  model  of  social 
interactions  between  independent  agents  and  utilizes  swarm  intelligence  to  achieve  the 
goal  of  optimization  problems.  It  was  chosen  here  since  the  associated  algorithm  has 
the  same  or  better  performances  comparing  to  other  global  optimization  programs. 

4.  Extension  of  the  program  to  include  the  analysis  of  spherical  cloaks.  By  this  it  will  be 
possible  to  analyze  double-curved  cloaking  structures. 
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The  main  realized  outcomes  of  the  project  are: 

•  We  have  developed  the  program  “UniaxCloak”  that  analyzes  cylindrical  and  spherical 
cloaks  made  from  the  uniaxial  materials.  The  program  can  fully  characterize  uniaxial 
cloaks,  i.e.  it  calculates  scattered  width  and  total  scattered  width  in  a  frequency  range 
of  interest.  In  the  cylindrical  case,  the  incident  wave  can  have  arbitrary  angle  of 
incidence. 

•  We  have  merged  the  program  “UniaxCloak”  with  the  global  optimization  program 
based  on  the  particle  swarm  optimization  (PSO)  method.  By  this  one  can  determine 
the  optimum  cloak  design  since  the  analysis  methods  based  on  the  transformation 
electromagnetics  do  not  take  all  the  effects  into  account. 

•  We  have  made  a  detailed  investigation  of  properties  of  cloaks  that  are  made  from 
uniaxial  materials.  In  practice,  such  cloaks  can  be  made  from  metamaterials,  for 
example  from  split-ring  periodic  structures  or  periodic  wire  structures. 
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Analysis  of  uniaxial  multilayer  cylinders 
used  for  invisible  cloak  realization 


3.1.  INTRODUCTION 

The  possibility  of  cloaking  objects  using  a  metamaterial  cover  has  extensively  been  studied 
(see  e.g.  [9]-[13]).  In  the  metamaterial  cloak  approach,  material  has  been  used  to  render  a 
volume  effectively  invisible  to  incident  radiation,  i.e.  to  squeeze  space  from  a  volume  into  a 
shell  surrounding  the  concealment  volume.  Coordinate  transformations  that  are  used  for  cloak 
design  do  not  influence  the  form  of  Maxwell’s  equations,  but  they  affect  permittivity  and 
permeability  tensors  (s  and  p,  respectively),  making  the  needed  materials  spatially  varying 
and  anisotropic.  When  viewed  externally,  the  concealed  volume  and  the  cloak  both  appear  to 
have  the  propagation  properties  of  free  space,  i.e.  they  appear  invisible  to  electromagnetic 
waves. 


Incident  wave 

^v\A — > 


Figure  3.1.  A  sketch  of  the  analyzed  structure. 

The  considered  structure  consists  of  a  PEC  cylinder  (object  being  cloaked)  and  of  a 
metamaterial  cloak.  A  sketch  of  the  analyzed  structure  is  given  in  Fig.  3.1.  For  the  cloak 
design,  a  coordinate  transformation  which  compresses  free  space  from  the  cylindrical  region  0 
<  r  <  b  into  the  concentric  cylindrical  shell  a  <  r  ’  <  b  is  applied,  where  a  and  b  represent  the 
cloak  inner  and  outer  radius,  respectively  (the  constitutive  structure  parameters  are  considered 
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in  the  cylindrical  coordinate  system).  The  used  transformation  leads  to  the  following 
expressions  for  the  components  of  the  permittivity  and  permeability  tensors  [9]: 


r-a 

£r=Mr= - 

r 


r 


(3.1a) 


Note  that  all  components  of  the  tensors  are  functions  of  radius,  which  implies  a  very 
complicated  metamaterial  design.  Such  a  cloak  is  referred  to  as  an  ideal  cloak  that  has  not  yet 
been  realized.  However,  several  cloak  designs  have  been  reported  which  claim  to  work 
properly  when  illuminated  with  a  normal  incident  plane  wave  of  specific  polarization  [14, 
15],  that  is,  with  either  the  electric  field  or  the  magnetic  field  parallel  to  the  z-axis  (TMZ 
polarization  or  TEZ  polarization,  respectively). 

The  metamaterial  cloak,  designed  by  Schurig  et  al.  [14],  is  intended  for  use  with  TMZ 
polarization  of  the  incident  wave  at  microwave  frequencies.  It  is  clear  that  only  ez,  //,  and  /v(/) 
components  are  relevant  in  this  case.  A  significantly  simplified  design  can  be  achieved  by 
fixing  values  of  sz  and  /iip  and  letting  only  /i,  to  vary  along  the  radial  direction: 


A,  = 1  • 


The  metamaterial  reported  in  [14]  used  the  unit  cell  of  approximately  cubic  shape,  filled  with 
one  split-ring  resonator.  The  realized  cloak  consisted  of  ten  concentric  cylinders,  each  of 
which  was  three  unit  cells  tall. 

In  the  simulations  presented  in  [14],  both  realizations  of  the  cloak  were  considered  (i.e., 
the  cloak  based  on  ideal  material  as  well  as  the  cloak  based  on  the  simplified  material  with 
spatially  varying  permeability).  When  simulating  the  simplified  version,  the  continuous 
variation  of  jur(r)  through  the  structure  was  approximated  by  a  10-step  piecewise  constant 
function,  that  mimicked  the  concentric  rings  of  the  fabricated  cloak.  The  results  revealed  that 
the  cloak  reduced  both  the  back-scattered  field  (reflection)  and  the  forward- scattered  field 
(shadow).  The  described  cloak  was  also  prototyped  and  embedded  into  scattering  chamber. 
The  distribution  of  the  vertical  component  of  the  electric  field  around  the  cloak  was  scanned, 
and  the  obtained  results  showed  some  decrease  of  the  scattered  field.  Although  this 
experiment  proved  the  basic  idea  of  cloaking,  the  level  of  scattered  field  has  not  been 
quantified. 

Similarly  to  the  realization  of  a  so-called  TMZ  cloak,  there  is  another  realization  of 
metamaterial  cloak  by  Cai  et  al.  [15],  This  cloak  is  intended  to  work  with  TEZ  polarized  wave 
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at  optical  frequencies.  In  TEZ  cloak  only  components  fiz,  sr  and  s,p  are  relevant,  thus  it  is  again 
possible  to  simplify  the  metamaterial  design  by  allowing  only  er  to  vary  in  radial  direction: 

Az=1 

/  \2  /  7  \2 


The  required  distribution  of  sr  is  realized  using  metal  wires  embedded  in  a  dielectric  material, 
i.e.  with  a  thin  wire  plasma-like  metamaterial.  Here,  the  effective  permittivity  varies  in  radial 
direction  and  its  real  part  has  to  exhibit  the  required  behaviour  (0<Re{  sr }  <  1 )  with  negligible 
imaginary  part  (i.e  with  negligible  losses)  [15]. 

The  purpose  of  this  chapter  is  to  provide  the  analysis  of  uniaxial  multilayer  cylinders  used 
for  invisible  cloak  realization.  Furthermore,  we  will  discuss  how  dispersion  influences  cloak 
properties,  i.e.  we  will  try  to  determine  the  frequency  bandwidth  of  the  considered  cloaks. 


r-a 


K  r  ) 


D 

\b-a  j 


\2 


(3.1c) 


3.2.  THEORETICAL  CONSIDERATIONS 

We  consider  a  circular-cylindrical  cloak  with  the  axis  in  z-direction  (see  Fig.3.1).  The 
incident  plane  wave  propagates  in  the  positive  x-direction  (i.e.  the  normal  incidence,  the 
oblique  incidence  will  be  treated  in  section  5),  while  the  electric  field  is  assumed  to  be 
parallel  with  the  z-axis  (TMZ  polarization;  TEZ  polarization  is  going  to  be  treated  afterwards). 
In  cylindrical  coordinates  the  E-  electric  field  component  is  given  in  the  form: 

oo 

£-(r,P)=  Z  j"L(Kr)e-J’r,  (3.2) 


where  Jm  is  the  m-th  order  Bessel  function  of  the  first  kind  and  ko  is  the  wave  number 
(k(,=2n/Xo).  Since  the  incident  wave  propagates  perpendicularly  to  the  cylinder  axis,  there  is 
no  variation  of  the  field  in  the  z-direction  (i.e.  d  / dz  =  0 ).  Consequently,  the  curl  Maxwell 
equations  can  be  written  as: 


1  dE__ 
jcof^/  d(p 


1  dE z 

dr 


]_ 

r 


aHr  ' 

d<p  , 


JojsEz. 


(3.3) 
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By  combining  these  equations  and  assuming  piecewise  constant  approximation  of  radial 
permeability  variation,  we  obtain  the  following  form  of  Bessel’s  differential  equation  for  the 
Ez  component  of  the  electric  field: 


2  d2Ez 
r  - r -  +  r 

dr 2 


^  +  (tV-m2)£z=0, 

dr 


(3.4) 


where  k  =  co2 ju(psz  and  m  =  m 


The  general  solution  of  (3.4)  is: 

c.j.cfco+cyOfe-),  (3.5) 

where  C i  and  C2  represent  constants  (to  be  determined  from  the  boundary  conditions),  while 
J-  and  H^2)  are  Bessel  functions  of  the  first  kind  and  Hankel  functions  of  the  second  kind  of 
order  m ,  respectively.  For  the  TEZ  cloaking  cylinder  we  obtain  the  equivalent  differential 
equation,  by  applying  the  substitutions  Ez-^>  Hz,  s— »p  and  p— Note  that  the  resulting  order 
of  Bessel  functions  ( m )  for  both  cloaks  is  not  integer  and  that  it  depends  on  both 
permeability  and  permittivity  tensors  (p  and  a). 

The  multilayered  structure  is  analyzed  using  G1DMULT  algorithm  [16],  which  calculates 
the  spectral-domain  Green's  functions  in  the  same  way  for  planar,  circular-cylindrical  and 
spherical  geometries.  The  algorithm  is  based  on  dividing  the  multilayer  problem  into 
equivalent  sub-problems,  one  for  each  dielectric  layer.  For  normal  incidence  the  fields  inside 
each  layer  of  the  cylindrical  structure  are  of  the  form  (e.g.  Ez  component): 

Ez(r,m)  =  almJm(klr)  +  blmll(2)(klr)  (3.6) 

Here  k‘  is  the  wave  number  in  z-th  layer,  and  a‘m  and  blm  are  the  unknown  coefficients  to  be 

determined.  The  equivalent  sub-problems  are  interrelated  (and  the  coefficients  alm  and  b'm 

are  calculated)  by  enforcing  the  continuity  of  tangential  components  of  electric  and  magnetic 
fields  over  the  layer  interfaces,  in  spectral  domain.  The  implementation  of  radially  anisotropic 
structures  into  the  G1DMULT  algorithm  is  quite  simple  -  one  only  needs  to  change  the  order 
of  Bessel/Hankel  function  in  each  anisotropic  layer  according  to  the  equation  (3.5).  More 
details  about  G1DMULT  can  be  found  in  [16]. 


3.3.  RESULTS  OF  THE  CLOAK  ANALYSIS 

In  order  to  characterize  the  level  of  “invisibility”  achieved  by  metamaterials,  we  have 
calculated  the  total  scattering  width  (gt)  of  the  cloaked  cylinder,  i.e.  the  total  scattered  power 
per  unit  length  normalized  with  the  incident  power  density  [17].  In  scientific  literature  there  is 
another  parameter  for  characterizing  scattering  objects  -  the  equivalent  blockage  width  (Weq). 
ft  is  a  complex-valued  parameter  (introduced  in  [6])  that  represents  the  width  of  an  ideal 
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shadow  which  produces  the  same  forward-scattered  field  as  the  cylinder  that  is  being 
observed  (in  our  case,  a  cloaking  cylinder).  Both  Weq  and  ct  are  quantities  that  actually  show 
how  wide  the  cylinder  seems  to  electromagnetic  waves.  There  is  a  simple  relation  between 
these  two  parameters  for  lossless  scatterers  -  ct  is  just  equal  to  2Re(Weq).  Therefore,  the 
considered  cloaks  can  be  characterized  with  total  scattering  width  only.  We  have  also 
calculated  the  angular  variation  of  bistatic  scattering  width  (020)  in  order  to  check  if  there  is 
any  direction  with  distinctively  stronger  scattered  field. 


3.3.1.  Ideal  cloak 

The  metamaterial  constitutive  parameters  have  been  calculated  by  equation  (3.1a).  Our 
intention  was  to  thoroughly  investigate  cloak  that  was  experimentally  realized  by  Schurig  et 
al.  [14].  Therefore,  we  have  selected  the  same  dimensions  and  the  working  frequency  as  in 
the  experimental  model  case.  The  cylinder  axis  is  set  in  z-direction  and  the  cloak  inner  and 
outer  radii  are  a  =  2.71  cm  and  b  =  5.89  cm,  respectively.  All  components  (r,  cp  and  z)  of 
permittivity  and  permeability  tensors  have  gradients  as  a  function  of  radius.  Such  anisotropy 
has  been  approximated  by  1-  to  10-  steps  piecewise  constant  functions  representing  the 
stepwise  cloak  realization,  in  order  to  analyze  how  the  subtlety  of  the  approximation  of  radial 
anisotropy  influences  the  total  scattering  width,  i.e.  “the  invisibility”. 

y  a 


& 


> 

X 


-aa/v — * 

Incident  plane  wave 


la 


Figure  3.2.  A  sketch  of  the  analyzed  structure.  The  cloak  contains  N  concentric  cylindrical 

metamaterial  layers. 

In  order  to  compare  different  cloak  realizations,  first  we  performed  one  referent  simulation 
for  the  case  without  any  cloak,  i.e.  with  only  PEC  circular  cylinder  present.  The  radius  of  the 
considered  cylinder  is  a  =  2.71  cm.  The  total  scattering  width  of  the  PEC  cylinder  is  equal  to 
12.71  cm  for  TMZ  and  8.98  cm  for  TEZ  polarization  of  the  incident  wave  at  the  working 
frequency  of  8.5  GHz,  see  Fig.  3.3.  By  calculating  the  total  scattering  widths  for  the  cases 
with  cloaks  present  it  also  becomes  possible  to  compare  the  obtained  invisibility  for  the  two 
analyzed  cases  (the  case  with  the  cloak  present  and  the  case  without  the  cloak  -  a  referent 
case).  The  obtained  invisibility  could  also  be  described  by  total  scattering  width  reduction 
(ratio  of  total  scattering  widths  for  the  two  cases). 
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Figure  3.3.  Total  scattering  width  of  PEC  cylinder  for  normal  incidence.  Both  TMZ  and  TEZ 

polarizations  are  shown. 


Due  to  the  full  anisotropy,  the  ideal  cloak  is  supposed  to  work  for  arbitrary  polarization. 
Therefore,  complete  analysis  is  provided  for  both  the  cloak  excited  by  the  TEZ  and  the  cloak 
excited  by  TMZ  polarized  waves.  The  simulations  have  been  performed  at  the  central 
frequency  (8.5  GHz).  The  calculated  total  scattering  width  and  the  angular  variation  of 
bistatic  scattering  width  are  shown  in  Figs.  3.4  -  3.7.  It  can  be  seen  that  by  increasing  the 
number  of  metamaterial  layers,  the  achieved  invisibility  improves.  Thus,  it  can  be  concluded 
that  the  10-layer  approximation  of  the  continuous  anisotropic  structure  is  good  enough  for 
practical  purposes  (even  5-layers  realization  has  a  very  good  performance).  It  can  also  be 
noted  that  the  level  of  scattered  field  is  low  for  both  polarizations  (Figures  3.6  and  3.7)  and, 
as  expected,  the  invisibility  was  achieved  for  both  excitations. 
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Figure  3.4.  Normalized  total  scattering  width  vs.  number  of  layers  for  ideal  cloak 

(TMZ  polarization). 


Figure  3.5.  Normalized  total  scattering  width  vs.  number  of  layers  for  ideal  cloak  (TEZ 

polarization). 
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Figure  3.6.  Electric-field  distribution  in  the  vicinity  of  the  10-layers  realization  of  ideal 

cloaking  cylinder  (TMZ  polarization). 


Figure  3.7.  Normalized  bistatic  scattering  width  of  10-layers  realization  of  ideal  cloaking 

cylinder  (TEZ  and  TMZ  polarization). 
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3.3.2.  TMZ  cloak  (Schurig  cloak) 

For  the  TMZ  cloak,  the  coordinate  transformations  for  metamaterial  design  have  been 
simplified,  so  the  metamaterial  constitutive  parameters  were  calculated  by  equation  (3.1b). 
Such  simplified  cloak  should  work  only  for  TMZ  polarization  and  thus  only  sz,  ji,  and 
components  are  relevant.  From  equation  (3.1b)  it  can  also  be  seen  that  both  the  permittivity 
tensor  component  in  z-direction  and  the  permeability  tensor  component  in  ^-direction  are 
constant.  Furthermore,  only  radial  anisotropy  of  permeability  tensor  (jur(r))  is  present.  This 
anisotropy  is  approximated  with  1  to  10-  steps  piecewise  constant  functions  in  a  way 
mentioned  before. 

The  analysis  of  how  the  subtlety  of  the  approximation  of  radial  anisotropy  influences  the 
total  scattering  width  (i.e.  the  “invisibility”)  has  been  given  only  for  excitation  by  TMZ 
polarization  (Fig.  3.8).  It  can  be  seen  that  for  structures  with  more  than  5  layers  there  is 
practically  no  improvement  of  the  obtained  invisibility  -  the  obtained  total  scattering  width 
reduction  is  around  3.  The  main  reason  for  such  a  small  gain  was  the  reflection  of  the  incident 
wave  from  the  cloak  surface  due  to  the  impedance  mismatch.  The  electric-field  distribution  in 
the  vicinity  of  the  TMZ  cloak  is  given  in  Fig.  3.9,  and  it  can  be  seen  that  the  reflection  from 
the  cloak  surface  causes  ripples  in  the  electric-field  distribution. 

For  the  cross-polarization  excitation  (TEZ)  the  radial  anisotropy  of  /u,  does  not  affect  the 
value  of  calculated  total  scattering  width.  This  is  because  the  magnetic  field  is  parallel  to  z- 
axis  and  thus  it  is  zero  in  radial  direction.  Fig.  3.10  shows  the  angular  variation  of  bistatic 
scattering  width  of  the  TMZ  cloak  for  both  polarizations.  Note  that  the  level  of  scattered  field 
for  the  TEZ  polarization  is  much  larger  than  the  one  calculated  for  the  TMZ  polarization. 
Therefore  the  TMZ  cloak  is  indeed  unsuitable  for  cross-polarization  excitation,  as  has  been 
presumed. 


Figure  3.8.  Normalized  total  scattering  width  vs.  number  of  layers  for  TMZ  cloak 

(TMZ  polarization). 
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Figure  3.9.  Electric-field  distribution  in  the  vicinity  of  the  10-layers  realization  of  the  TMZ 

cloaking  cylinder  (TMZ  polarization). 
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Figure  3.10.  Normalized  bistatic  scattering  width  of  10-layers  realization  of  TMZ  cloaking 

cylinder  (TMZ  and  TEZ  polarization). 
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3.3.3.  TEZ  cloak  (Cai  cloak) 

The  TEZ  cloak  was  designed  as  a  completely  dual  cloak  to  the  TMZ  cloak  described  in  the 
previous  section  and  the  metamaterial  constitutive  parameters  were  calculated  by  equation 
(3.1c).  Such  simplified  cloak  is  thus  designed  to  work  only  for  TEZ  polarization  and  only  juz,  er 
and  c(p  components  of  the  associated  tensors  are  relevant.  Complementary  to  the  TMZ  cloak, 
for  the  TEZ  cloak  only  radial  variation  of  the  electric  permittivity  ( £,(>'))  is  present.  This 
anisotropy  is  also  approximated  with  1  to  10-  steps  piecewise  constant  functions. 

The  results  of  the  simulations  performed  for  the  TEZ  cloak  are  shown  on  Figs.  3.1 1  and  3.12. 
As  expected,  the  results  are  complementary  to  the  ones  provided  for  TMZ  cloak.  The  realized 
total  scattering  width  reduction  is  again  smaller  than  in  the  ideal  case,  i.e.  it  is  around  3  for 
structures  with  more  than  5  layers.  As  expected,  the  cloak  is  found  to  be  suitable  only  for  TEZ 
excitation. 


Figure  3.11.  Normalized  total  scattering  width  vs.  number  of  layers  for  TEZ  cloak 

(TEZ  polarization). 
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Figure  3.12.  Normalized  bistatic  scattering  width  of  10-layers  realization  of  TEZ  cloaking 

cylinder  (TMZ  and  TEZ  polarization). 


3.4.  LIMITATIONS  OF  SIMPLIFIED  CLOAK  REALIZATIONS 


The  performed  simulations  show  that  for  the  simplified  cloak  realizations  (TMZ  and  TEZ 
cloaks)  the  achieved  total  scattering  width  reduction  is  around  3,  though  only  for  one  certain 
excitation.  In  reality,  the  magnetic  permeability  and  the  electric  permittivity  in  metamaterials 
are  always  frequency  dependent.  Therefore,  the  frequency  variations  of  //  and  e  would  have  an 
effect  on  the  level  of  achieved  invisibility  (see  eq.  3.4  and  3.5). 

The  frequency  dependence  of  the  magnetic  permeability  (realized  with  some  kind  of  split 
rings)  is  given  by  the  so-called  Lorentz  model  [18]: 


/k//=|- 


f  2-f2 

J  mp  J  0 
/2-/o2-J>/' 


(3.7) 


Here /is  the  frequency  of  the  signal,/^  denotes  the  frequency  at  which  /ief/= 0  for  the  lossless 
case  (null-point  of  the  function),  fo  is  the  frequency  at  which  jiefy  diverges  (the  pole  of  the 
function),  while  the  factor  y  represents  the  losses.  In  our  calculations  we  have  approximated 
fmP  by  fmp  =  1.02  fo,  following  the  measured  values  of  practical  SRR-based  metamaterials 
published  in  [18]  and  [19]  (see  Figure  3.13). 
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Figure  3.13.  Frequency  dependence  of  /ur  (from  [18]  and  [19]). 


On  the  other  hand,  the  frequency  dependence  of  the  electric  permittivity  (realized  with 
some  kind  of  array  of  thin  wires)  is  given  by  the  so-called  Drude  model  [18]: 


£eff  ~ 


1  — 


f-jyf 


'  host  * 


(3.8) 


Here  /  is  the  frequency  of  the  signal  and  fp  represents  the  frequency  at  which  seff  =  0  for  the 
lossless  case  (null-point  of  the  function).  Factor  y  represents  the  losses,  while  Shost  is  the 
permittivity  of  the  medium  which  is  hosting  the  metamaterial  structure. 

Frequency  dependencies  of  relative  magnetic  permeability  and  relative  electric  permittivity 
in  radial  direction  are  shown  in  Fig.  3.9.  For  simplicity,  only  the  lossless  case  is  being 
considered  here  (i ,e.jyf=  0). 

For  the  central  frequency  /  =  8.5  GHz  it  is  assumed  that  the  radial  components  of  the 
permittivity  and  permeability  tensors  are  jur=  0.14  (for  TMZ  cloak)  and  er  =  0.48  (for  TEZ 
cloak).  These  are  the  values  needed  for  the  cloak  realizations  with  one  layer  only.  It  is  now 
obvious  that  if  //,•  or  sr  throughout  the  structure  are  designed  according  to  eq.  3.1b  and  3.1c, 
these  values  can  be  achieved  only  at  some  certain  frequency.  Due  to  dispersion,  the  values  of 
Hr  or  sr  will  grow  with  increasing  frequency.  Therefore,  the  achieved  level  of  invisibility  will 
also  change. 
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Figure  3.14.  Frequency  dependence  of  / ir  and  er  (Lorentz  and  Drude  model). 

The  simulations  considered  the  simplified  cloak  realizations,  in  order  to  calculate  the  total 
scattering  width  with  dispersion  included.  The  cloaks  were  designed  with  ten  layers  of 
metamaterial,  which  is  a  good  approximation  of  continuous  radial  change  of  permittivity  or 
permeability  values,  as  shown  in  section  3.2.  The  analysis  is  provided  in  a  narrow  frequency 
range  about  the  central  frequency  (8.5  GHz),  for  which  the  values  of  /j,r  or  sr  had  been 
designed.  Figs.  3.15  and  3.16  show  the  frequency  dependence  of  the  invisibility  parameters 
for  the  simplified  cloaks  -  the  TMZ  and  TEZ  cloak,  respectively.  The  comparison  with  the 
referent  case  (PEC  cylinder  without  any  cloak)  is  also  given. 

We  can  define  the  bandwidth  of  the  cloak  as  the  frequency  range  (normalized  to  the  central 
frequency),  in  which  the  total  scattering  width  is  smaller  than  the  total  scattering  width  of  the 
hidden  object  (PEC  cylinder  in  our  case).  It  can  be  seen  that  the  invisibility  has  been  achieved 
only  in  a  narrow  frequency  range  around  the  central  frequency  (8.5  GHz).  Therefore,  it  has 
also  been  found  that  the  metamaterial  cloak  is  not  suitable  for  applications  that  require  a 
larger  frequency  bandwidth  (e.g.  invisibility  to  radar  signal).  If  we  suppose  that  the  cloak 
behaves  similarly  at  frequencies  lower  than  the  central  one,  we  find  that  the  approximate 
bandwidths  of  the  TMZ  and  TEZ  cloaks  are  0.02  GHz  and  0.2  GHz,  respectively  (0.24%  and 
2.4  %  of  the  central  frequency).  The  TEZ  cloak  has  a  larger  bandwidth  due  to  the  weaker 
frequency  dependence  of  constitutive  parameters  (see  Fig.  3.14). 
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Figure  3.15.  Normalized  total  scattering  width  vs.  frequency  for  TMZ  cloak  with  included 

dispersion  (TMZ  polarization). 


Figure  3.16.  Normalized  total  scattering  width  vs.  frequency  for  TEZ  cloak  with  included 

dispersion  (TEZ  polarization). 


The  bandwidth  results  of  the  TMZ  cloak  strongly  depend  on  the  ratio  between  the  frequencies 
fmp  (the  frequency  at  which  juefj=0  for  the  lossless  case)  and  fo  (the  frequency  at  which  /jejf 
diverges  for  the  lossless  case).  Since  this  ratio  depends  on  the  particular  realization  of  the 
split-ring  resonator,  we  have  investigated  the  dependence  of  the  cloak  bandwidth  on  that 
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ratio.  In  Figure  3.17  we  have  illustrated  this  dependence.  It  can  be  seen  that  with  the  ratio 
=  1.5  fo  one  can  obtain  10  times  large  bandwidth,  i.e.  bandwidth  of  around  2.4  %  for  the  TMZ 
polarized  incident  wave.  However,  larger  ratio  in  practice  means  larger  losses  and  worse 
cloak  behavior  (the  dependence  of  the  cloak  on  losses  will  be  discussed  at  the  end  of  this 
section).  Therefore,  one  can  conclude  that  if  the  cloak  is  made  from  resonant  type  of  elements 
the  bandwidth  will  be  quite  narrow. 


Frequency  (GHz) 


Figure  3.17.  Normalized  total  scattering  width  as  a  function  of  ratio  /m/;  /  f,  for  TMZ  cloak 

(TMZ  polarization). 

In  order  to  verify  our  bandwidth  investigation  (and  to  verify  the  developed  code)  we  have 
compared  the  results  calculated  with  the  developed  program  and  the  results  found  in  scientific 
literature.  In  [20]  the  bandwidth  of  the  cloak  made  from  elements  with  the  Lorentz  dispersion 
behavior  was  investigated.  The  authors  have  analyzed  the  10-layer  Schurig  cloak,  and  they 
have  supposed  the  non-uniform  distribution  of  the  fmp  /  fo,  see  table  3.1.  The  obtained 
bandwidth  is  around  2  %.  The  comparison  of  the  total  scattering  width,  calculated  with  our 
program  and  with  the  finite-element  solver  Comsol  Multiphysics,  is  given  in  Fig.  3.18.  The 
agreement  with  the  results  is  very  good. 

The  influence  of  losses  of  the  periodic  structure  is  investigated  in  Fig.  3.19  where  the  total 
scattering  width  is  calculated  as  a  function  of  losses  (described  with  factor  y  from  equation 
(3.7)).  It  can  be  seen  that  if  the  value  of  y  is  smaller  than  0.001  fo  then  the  difference  between 
lossy  and  lossless  cases  is  negligible.  However,  for  larger  losses  the  cloak  loose  its  cloaking 
behavior  and  it  starts  to  behave  as  an  absorber  (i.e.  the  forward  scattered  field  is  drastically 
enlarged  and  the  backward  scattered  field  is  reduced).  That  can  be  clearly  seen  in  Figure  3.20 
which  illustrate  the  electric  field  distribution  in  the  vicinity  of  the  10-layers  realization  of  the 
TMZ  cloaking  cylinder  with  y  =  0.01  fo  .  One  can  clearly  see  that  the  cloak  starts  to  behave  as 
a  absorber  resulting  in  strong  forward  scattered  field. 
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Shell 

blip 

fo 

fmp/f0 

1 

8.4901 

3.98772 

2.1291 

2 

8.4344 

4.89096 

1.7245 

3 

8.3657 

5.35577 

1.5620 

4 

8.2953 

5.63920 

1.4710 

5 

8.2237 

5.80921 

1.4156 

6 

8.1597 

5.95249 

1.3708 

7 

8.0932 

6.02870 

1.3424 

8 

8.0336 

6.09923 

1.3172 

9 

7.9804 

6.16397 

1.2947 

10 

7.9258 

6.20024 

1.2783 

Table  3.1.  The  dispersion  parameters  of  the  cloak  investigated  in  [20], 


50 

45 

_40 

E 

£-35 

|  30 
o> 

•c  25 
ju 

S  20 

in 

1 15 

H 

10 

5 

0 

8 


O  G1DMULT 
*  General  solver 

— i - 1 - 1 - 

-  $  . | . 

$ 

* 

. ! . I .  *  j 

. 1 . I . i . *  i  ?- 

. ! . 1 .  . 

$ 

9 

. f  # 

$  * 

♦ 

♦ 

*  ♦ 

a  _ _ i _ i _ 

8.4  8.5  8.6 

Frequency  (GHz) 


8.7 


8.8 


Figure  3.18.  Comparison  of  the  total  scattering  width  calculated  by  the  developed  program 
(based  on  the  G1DMULT  algorithm)  and  by  the  general  finite-element  solver  Comsol 
Multiphysics  [20],  The  dispersion  parameters  of  the  cloak  are  given  in  table  3.1. 
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Figure  3.19.  Electric-field  distribution  in  the  vicinity  of  the  10-layers  lossy  realization  of  the 
TMZ  cloaking  cylinder  (TMZ  polarization).  The  losses  are  y  =  0.01 -fo  . 
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4  OPTIMIZATION  OF  UNIAXIAL  MULTILAYER 
CYLINDERS  USED  FOR  INVISIBLE  CLOAK 

REALIZATION 
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Optimization  of  Uniaxial  Multilayer  Cylinders 
Used  for  Invisible  Cloak  Realization 


4.1  Global  optimization  techniques 

Global  optimization  is  a  branch  of  applied  mathematics  and  numerical  analysis  that  deals 
with  the  optimization  of  a  function  or  a  set  of  functions  to  some  criteria.  They  can  be 
formulated  as  a  TV-dimensional  minimization  problem  as  follows: 

min/(x),  x  =  [x1,x2,...,xAf ] ,  (4.1) 

where  N  is  the  number  of  the  parameters  to  be  optimized.  There  are  many  global  optimization 
techniques  available,  but  so  far  in  electromagnetics  only  few  have  extensively  been  used. 
They  are  the  genetic  algorithm,  the  particle  swarm  optimization  (PSO),  the  simulated 
annealing,  the  ant  colony  optimization  and  the  multidimensional  conjugate  gradient  method. 
All  the  global  optimization  techniques  prove  to  be  usable,  but  at  this  point  we  have 
concentrated  on  the  PSO.  We  have  chosen  it  since  the  associated  algorithm  is  rather  easy  and 
straightforward  to  implement,  the  parameters  used  to  tweak  its  performance  are  easily 
understandable,  and  it  has  the  same  or  better  performance  compared  to  other  global 
optimization  algorithms.  We  have  implemented  three  different  kinds  of  PSO  algorithm  in 
order  to  compare  their  performance  on  optimization  of  complex  electromagnetic  problems, 
i.e.  synthesis  of  various  multilayered  electromagnetic  cloak  designs. 

4.1.1.  Classical  particle  swarm  optimization 

The  classical  particle  swarm  optimization  (PSO)  is  a  robust  stochastic  evolutionary 
computation  technique  based  on  the  movement  and  intelligence  of  swarms.  It  was  initially 
introduced  by  J.  Kennedy  and  R.  C.  Eberhart  in  1995  [21],  but  in  electromagnetic  community 
it  has  not  been  extensively  used  until  2004  when  J.  Robinson  and  Y.  Rahmat-Samii  published 
their  paper  [22], 

As  authors  [22]  have  conveniently  put  it,  the  PSO  can  best  be  understood  through  an 
analogy  similar  to  the  one  that  led  to  the  development  of  the  PSO.  Imagine  a  swarm  of  bees  in 
a  field.  Their  goal  is  to  find  in  the  field  the  location  with  the  highest  density  of  flowers. 
Without  any  knowledge  of  the  field  a  priori,  the  bees  begin  in  random  locations  with  random 
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velocities  looking  for  flowers.  Each  bee  can  remember  the  locations  that  it  found  the  most 
flowers,  and  somehow  knows  the  locations  where  the  other  bees  found  an  abundance  of 
flowers.  Tom  between  returning  to  the  location  where  it  had  personally  found  the  most 
flowers  and  exploring  the  location  reported  by  others  to  have  the  most  flowers,  the  ambivalent 
bee  accelerates  in  both  directions  altering  its  trajectory  to  fly  somewhere  between  the  two 
points  depending  on  whether  nostalgia  or  social  influence  dominates  its  decision.  Along  the 
way,  a  bee  might  find  a  place  with  a  higher  concentration  of  flowers  than  it  had  found 
previously.  It  would  then  be  drawn  to  this  new  location  as  well  as  the  location  of  the  most 
flowers  found  by  the  whole  swarm.  Occasionally,  one  bee  may  fly  over  a  place  with  more 
flowers  than  had  been  encountered  by  any  bee  in  the  swarm.  The  whole  swarm  would  then  be 
drawn  toward  that  location  in  additional  to  their  own  personal  discovery.  In  this  way  the  bees 
explore  the  field:  overflying  locations  of  greatest  concentration,  then  being  pulled  back 
toward  them.  Constantly,  they  are  checking  the  territory  they  fly  over  against  previously 
encountered  locations  of  highest  concentration  hoping  to  find  the  absolute  highest 
concentration  of  flowers.  Eventually,  the  bees’  flight  leads  them  to  the  one  place  in  the  field 
with  the  highest  concentration  of  flowers.  Soon,  all  the  bees  swarm  around  this  point.  Unable 
to  find  any  points  of  higher  flower  concentration,  they  are  continually  drawn  back  to  the 
highest  flower  concentration.  In  attempting  to  model  this  behavior,  Kennedy  and  Eberhart 
realized  that  they  had  stumbled  upon  an  optimizer. 


& 
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Figure  4.1.  Bees  in  a  search  of  field  for  flowers.  They  are  attracted  to  location  with  highest 
concentration  of  flowers  they  have  found  and  to  the  location  with  the  highest  concentration  of 
flowers  that  the  entire  swarm  has  found  (a).  Eventually  all  the  bees  converge  to  the  region 
where  the  highest  concentration  of  flowers  that  the  entire  swarm  has  found  is  (b).  (image 
taken  from  [22]) 


The  language  used  to  discuss  the  PSO  follows  from  the  analogy  of  particles  in  a  swarm: 

(1)  Particle  or  Agent:  Each  individual  in  the  swarm  (bee)  is  referred  to  as  a  particle  or 
agent.  Particles  act  individually:  accelerating  toward  the  best  personal  and  best 
overall  location  while  checking  the  value  of  its  current  location  in  each  time-step. 

(2)  Position:  In  the  analogy  above  position  is  the  bee’s  place  in  the  field.  In  this  case, 
position  is  represented  by  coordinates  on  the  plane.  In  general  this  idea  can  be 
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extended  into  any  iV-dimensional  space  according  to  the  problem  at  hand.  This  N- 
dimensional  space  is  the  solution  space  for  the  problem  being  optimized,  where  any 
set  of  coordinates  (position  vector)  represents  a  solution  to  the  problem.  In  general 
these  can  be  any  values  needed  to  be  optimized.  Reducing  the  optimization  problem 
to  a  set  of  values  that  could  represent  a  position  in  solution  space  is  an  essential  step 
in  utilizing  the  PSO. 

(3)  Fitness:  As  in  all  evolutionary  computation  techniques  there  must  be  some  function 
or  method  to  evaluate  the  goodness  (fitness)  of  a  position.  The  fitness  function  takes 
the  position  in  the  solution  space  (position  vector)  and  returns  a  single  number 
representing  the  value  of  that  position.  In  the  bees  searching  for  flowers  analogy,  the 
fitness  function  would  simply  be  the  flower  density:  the  higher  the  density,  the  better 
the  location.  The  fitness  function  provides  the  interface  between  the  physical  problem 
and  the  optimization  algorithm.  Special  care  has  to  be  taken  into  account  if  multigoal 
analysis  is  carried  out.  More  on  this  will  be  said  in  subsequent  chapters. 

(4)  Personal  best  (pbest):  In  the  bees  analogy,  each  bee  remembers  the  location  where  it 
personally  encountered  the  most  flowers.  This  location  with  the  highest  fitness  value 
personally  discovered  by  a  bee  is  known  as  the  personal  best  or  pbest.  Each  bee  has 
its  own  personal  best  determined  by  the  path  that  it  has  flown.  At  each  time-step,  the 
bee  checks  if  the  current  location  has  a  higher  fitness  value  and  if  it  does,  the  personal 
best  is  replaced  with  bee’s  current  location. 

(5)  Global  best  (gbest):  Each  bee  also  knows  the  highest  concentration  of  flowers 
discovered  by  the  entire  swarm.  This  location  of  highest  fitness  encountered  is  known 
as  the  global  best  or  gbest.  For  the  entire  swarm  there  is  one  global  best  to  which  all 
bees  are  attracted.  At  each  point  along  their  path  bees  compare  the  fitness  of  their 
current  location  to  that  of  global  best.  If  any  bee  is  at  a  location  of  higher  fitness, 
global  best  is  replaced  by  that  bee’s  current  position. 


The  algorithm  goes  as  follows: 

(1)  Define  the  solution  space 

When  trying  to  optimize  physical  problems,  one  needs  to  select  the  parameters  of  the 
problem  that  are  suitable  for  optimization.  After  selecting  N  parameters,  for  each 
parameter  reasonable  (physical)  range  in  which  to  search  for  optimal  solution  must  be 
given.  For  each  dimension  in  an  TV-dimensional  optimization  Xmin  and  Xmax  values  have 
to  be  specified. 

(2)  Define  fitness  function 

The  fitness  function  is  the  link  between  the  optimization  algorithm  and  the  physical 
problem  that  is  to  be  optimized.  The  fitness  function  takes  N  input  values  (for  each 
dimension,  i.e.  parameter  to  be  optimized)  and  produces  one  single  value  that  should 
unambiguous  define  the  quality  of  the  selected  solution.  In  case  of  multi-goal  analysis, 
weighting  factors  are  typically  introduced  -  so  each  goal  can  be  weighted  differently 
depending  on  the  importance.  For  example,  if  we  are  trying  to  maximize  the  antenna 
gain  while  keeping  low  side  lobe  levels,  the  corresponding  fitness  function  could  be 
defined  as: 


Fitness  =  a  ■  Gain  -  j3  ■  SL(H )  -  y  ■  SL(E). 


(4.2) 
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Where  Gain  is  the  antenna  gain  at  broadside,  and  SL(H)  and  SL(E)  are  the  highest  side 
lobe  levels  in  H-  and  E-planes,  respectively,  a,  /?,  and  y  are  weighting  coefficients  that 
can  be  adjusted  to  maximize  the  gain,  or  produce  a  beam  with  low  side  lobes.  Their 
values  could,  for  instance,  be  1,  0.3  and  0.3,  and  by  this  way  the  optimizer  would  try 
to  maximize  the  gain,  but  keep  track  of  the  sidelobes. 

(3)  Define  initial  particle  locations  and  velocities 

Each  particle  begins  its  search  for  optimal  solution  at  a  random  location  with  a 
random  velocity  (direction  and  magnitude).  These  initial  values  are  preselected,  most 
often  by  random  selection,  but  can  also  be  defined  according  to  some  grid.  When 
particle  is  given  its  initial  location,  personal  best  of  each  particle  is  calculated  and 
updated,  and  among  them  the  initial  global  best  is  selected. 

(4)  Run  the  optimization  loop 

The  particles  are  systematically  flown  through  the  solution  space.  Each  particle  is 
moved.  The  movement  is  determined  by  the  velocity  vector.  Once  moved,  the  fitness 
is  calculated,  personal  and  global  best  are  updated  if  the  calculated  fitness  is  greater 
than  one  remembered  from  before  and  speed  is  updated.  The  equation  for  updating  the 
velocity  will  be  defined  later,  but  basically  the  velocity  is  determined  by  velocity  in 
previous  step  and  random  influence  of  the  distances  between  the  current  location  and 
personal  and  global  best  positions.  The  loop  is  repeated  until  some  stop  criterion  is 
met.  The  stop  criterions  can  be  many.  The  optimization  can  stop  once  certain  fitness 
level  is  reached,  if  predetermined  number  of  iterations  has  been  carried  out,  or  if  the 
optimization  has  not  improved  the  global  best  in  some  time.  To  more  easily  follow  the 
optimization  process  a  flow  chart  is  drawn. 


/ 
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Figure  4.2.  Flow  chart  of  the  PSO  algorithm,  (a)  The  whole  flow  chart.  The  part  inside 
dashed  rounded  rectangle  is  independent  of  the  physical  system  being  optimized,  (b)  The 
main  PSO  loop 
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UPDATE  VELOCITY 

The  velocity  of  the  particle  is  changed  according  to  its  position  relative  to  the  locations  of 
personal  and  global  best.  The  particle  is  accelerated  in  the  directions  of  these  locations  of 
greatest  fitness  according  to  the  following  equation: 

Vn  =  W' ■ V«  +  Cr  rand 0  •  (Pbest,n  ~  Xn  ) +  C2  '  ™nd 0  '  (gbes,,n  ~  Xn  )  »  (4.3) 

where  vn  is  the  velocity  of  the  particle  in  «-th  dimension  and  xn  is  the  particle  coordinate  in 
the  same  dimension.  The  vn  is  calculated  for  every  dimension  of  the  N  dimensional  solution 
space.  From  this  equation  it  can  be  seen  that  new  velocity  is  equal  to  old  velocity  scaled  by 
the  w  (called  the  inertia  factor)  and  increased  toward  the  locations  of  personal  and  global  best. 
ci  and  C2  are  scaling  factors  that  determine  the  relative  “pull”  of  personal  best  and  global  best. 
These  are  sometimes  referred  to  as  the  cognitive  and  social  rates,  respectively.  Increasing  c; 
encourages  exploration  of  the  solution  space  as  each  particle  moves  toward  its  own  personal 
best,  while  increasing  C2  encourages  exploitation  of  the  supposed  global  maximum.  The  three 
parameters  w,  c;  and  C2  can  be  used  to  fine  tune  the  search  process.  More  on  that  will  follow 
later.  The  random  number  function  rand(-)  returns  a  number  between  0.0  and  1.0.  This 
introduction  of  a  random  element  into  the  optimization  is  intended  to  simulate  the  slight 
unpredictable  component  of  natural  swarm  behavior,  together  with  initialization  process  that 
varies  the  starting  positions  and  velocities  of  the  particles. 

UPDATE  POSITION 

After  particle  fitness  has  been  evaluated  and  velocity  updated  accordingly,  it  is  time  to 
move  the  particle  to  a  new  location.  The  velocity  is  applied  for  a  given  time-step  At,  usually 
chosen  to  be  one  and  new  coordinate  is  computed  for  each  of  the  dimensions  according  the 
following  equation: 


x  =x, +At-v,  (4.4) 

SELECTION  OF  PARAMETERS 

Compared  to  some  other  global  optimization  techniques,  PSO  has  only  a  couple  of 
parameters  that  have  to  be  selected  prior  to  optimization.  The  goal  was  to  develop  an 
algorithm  with  an  optimal  balance  between  global  exploration  and  exploitation  of  local 
maxima,  i.e.  an  algorithm  that  would  quite  steadily  converge,  but  at  the  same  time  be  able  to 
escape  local  extremes.  The  first  problem  was  control  over  the  search  space.  Without 
boundaries  and  limits  on  the  velocities  and  due  to  the  random  nature,  many  particles  would 
fly  out  of  the  physically  meaningful  solution  space.  This  would  lead  to  loss  of  time,  as  the 
algorithm  would  still  calculate  the  fitness  (typically  computationally  most  expensive  part  of 
the  optimization  procedure)  of  these  particles  although  their  values  as  solutions  where 
practically  useless  if  the  search  space  was  defined  properly  prior  to  optimization  start. 

This  problem  was,  at  first,  tackled  by  introduction  of  velocity  clamping,  i.e.  enforcing  a 
maximum  allowed  velocity  called  Vmax.  Most  authors  agreed  that  it  is  best  to  set  Vmax  in  each 
dimension  to  the  dynamic  range  of  that  dimension.  Second  approach  includes  varying  value 
of  inertial  weight  in  an  attempt  to  strike  a  balance  between  global  exploration  and  local 
exploitation.  Larger  values  of  inertial  weight  tend  to  encourage  global  exploration  as  a  result 
of  the  particle  being  less  moved  by  the  pull  of  personal  and  global  best,  while  with  lower 
values  particles  are  more  rapidly  attracted  to  the  personal  and  global  best  positions.  Realizing 
the  importance  of  exploration  early  in  the  run,  and  the  increasing  importance  of  exploitation 
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of  maxima  as  the  run  progresses,  it  was  suggested  to  vary  the  inertial  weight  linearly  from  0.9 
to  0.4  over  the  course  of  the  optimization  run.  For  this  technique  to  be  successful,  it  is 
important  to  have  sufficient  number  of  iterations,  so  both  search  states  could  be  utilized 
sufficiently. 

Population  size  is  another  parameter  that  should  carefully  be  selected.  Large  populations, 
while  providing  the  most  thorough  exploration  of  the  solution  space,  come  at  the  cost  of  more 
fitness  evaluations  and  computation  time.  For  the  PSO,  previous  studies  claimed  that 
relatively  small  population  sizes  can  sufficiently  good  explore  a  solution  space  while  avoiding 
excessive  fitness  evaluations.  Our  results  show  that  small  population  is  not  necessarily  better 
choice  [24],  In  study  that  we  carried  out  while  optimizing  layered  circular-cylindrical 
dielectric  lens  antennas  we  found  out  that  sometimes  it  is  beneficial  to  increase  the  particle 
count  and  reduce  the  number  iterations.  In  this  way,  total  number  of  evaluations  is  kept  at 
same  levels,  but  better  optimization  results  can  be  obtained.  Increasing  the  particle  count  is 
especially  beneficiary  when  optimizing  complex  multidimensional  problems.  In  that  case, 
there  is  a  high  chance  that  the  particles  will  be  stuck  in  local  extremes.  Larger  number  of 
particles  that  communicate  to  each  other  are  more  likely  to  bypass  this  problem. 

BOUNDARY  CONDITIONS 

The  definition  of  the  solutions  space  is  one  of  the  first  steps  carried  out  while  preparing  a 
specific  problem  to  be  optimized  by  the  PSO.  If  the  solution  space  is  well  defined,  the 
optimization  procedure  will  evaluate  only  practically  realizable  candidates.  The  question  that 
naturally  arises  is:  what  happens  when  particle,  by  the  randomness  nature  of  its  velocity, 
comes  to  the  border  of  the  solution  space.  To  address  this  problem,  the  inventors  of  the  PSO 
have  imposed  three  different  boundary  conditions:  absorbing,  reflecting  and  invisible  walls. 

(1)  Absorbing  Walls:  When  a  particle  hits  the  boundary  in  one  of  the  dimensions,  the 
velocity  in  that  dimension  is  set  to  zero.  The  particle  does  not  stay  at  the  boundary  for 
long,  as  it  is  pulled  back  toward  the  allowed  solution  space  via  random  nature  of  the 
velocity  vector.  We  can  say  that  boundary  walls  absorb  the  energy  of  particles  trying 
to  escape  the  solution  space. 

(2)  Reflecting  Walls:  When  a  particle  hits  the  boundary  in  one  of  the  dimensions,  the 
sign  of  the  velocity  in  that  dimension  is  changed,  i.e.  the  particle  is  reflected  back 
toward  the  solution  space  with  no  loss  of  “kinetic”  energy. 

(3)  Invisible  Walls:  The  particles  are  allowed  to  fly  anywhere  without  any  physical 
restriction,  but  for  particles  that  roam  outside  the  allowed  solution  space  the  fitness  is 
not  calculated.  For  nearly  all  engineering  applications,  the  computationally  more 
expensive  portion  of  the  algorithm  is  the  fitness  evaluation,  so  by  using  the  invisible 
walls  the  execution  time  is  not  increased  significantly.  The  motivation  behind  this 
technique  is  to  save  computation  time  by  only  evaluating  what  is  in  the  allowed 
solution  space,  while  not  interfering  with  the  natural  motion  of  the  swarm. 

Most  authors  agree  that  the  “invisible  walls”  technique  provides  same  or  slightly  better  results 
than  other  wall  techniques.  For  this  reason  we  have  implemented  the  “invisible  walls” 
technique  in  our  optimization  routines. 
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absorbing  walls  reflecting  walls  invisible  walls 

Figure  4.3.  Graphical  interpretation  of  various  boundary  conditions  used  in  PSO. 


4.1.2.  Particle  swarm  optimization  with  local  best  topology 

In  the  classical  (traditional)  particle  swarm  described  above,  each  individual  has  insight  to 
all  particles  present  in  the  swarm,  i.e.  it  knows  and  uses  the  global  best  solution  to  determine 
its  velocity  vector  and  consequently  the  trajectory  by  which  it  searches  the  solution  space. 
This  represents  an  oversimplification  of  the  social-psychological  view  that  individuals  are 
more  affected  by  those  who  are  more  successful,  persuasive,  or  otherwise  prestigious.  In  the 
human  society,  it  is  more  accurate  to  say  that  the  social  neighborhood  provides  a  wealth  of 
possible  models  whose  behavior  may  be  emulated,  and  individuals  seem  to  be  affected  by 
some  kind  of  statistical  summary  of  the  state  of  their  immediate  social  network  rather  than  the 
unique  performance  of  one  individual  [25],  This  is  the  reason  for  introducing  different 
neighborhood  topologies. 

The  classical  PSO  algorithm  outlined  in  previous  section  used  a  kind  of  topology  that  is 
known  under  name  global  best  (or  gbest).  Each  particle,  in  this  topology,  is  influenced  by  the 
best  performing  individual  in  the  entire  population.  This  is  equivalent  to  a  social  network 
where  every  particle  is  connected  to  every  other  particle.  This  kind  of  topology  was 
acceptable,  at  least  according  to  Kennedy  and  Mendes  [25],  in  first  applications  where 
function  landscape  is  largely  made  up  of  long  gradients,  where  problem  is,  first,  to  find  the 
best  gradient  region  of  the  search  space  and,  second,  to  find  the  extreme  of  that  region. 

Many  problems,  including  the  one  studied  here  (synthesis  of  various  multilayered 
electromagnetic  cloak  designs),  contain  cliffs,  variable  interactions  and  other  features  that  are 
not  typified  by  smooth  gradients  -  for  this  kind  of  problems,  a  more  robust  algorithm  is 
needed,  at  least  according  to  [25],  A  new  topology,  termed  local  best  (or  lbest)  was  proposed 
in  order  to  deal  with  these,  more  difficult,  problems. 

In  local  best  topology,  the  subpopulations  can  search  diverse  regions  of  the  problem  space. 
One  part  of  the  population  can  concentrate  on  one  local  optimum,  while  other  part  can  search 
around  different  local  optimum.  The  flow  of  information  between  particles  is  slowed  down, 
and  that  can  prevent  premature  convergence  from  happening.  Further  tweaking  of  the  search 
can  be  made  by  varying  the  size  of  local  neighborhood,  i.e.  the  number  of  particles  that  each 
particle  communicates  to. 

To  summarize,  the  global  best  topology  (i.e.,  the  biggest  neighborhood  possible)  often 
converges  more  quickly  compared  to  the  local  best  topology,  but  at  the  same  time  it  is  also 
more  susceptible  to  the  attraction  of  local  optima  since  the  population  rushes  unanimously 
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toward  the  first  good  solution  found.  In  our  optimizations,  we  have  typically  set  the  size  of 
the  swarm  at  100,  following  the  reasoning  brought  forward  in  previous  section,  while  the 
local  neighborhood  was  set  at  5  and  was  defined  during  the  first  evaluation  based  on 
normalized  distance  between  particles. 


Figure  4.4.  Different  types  of  topologies,  (a)  comparison  between  gbest  topology  used  in 
classical  PSO  (every  particle  connects  to  all  others)  and  one  version  of  lbest  topology  where 
each  particle  is  connected  to  two  neighbors  (b)  von  Neumann  or  square  topology  -  each 
individual  is  connected  to  four  other  particles  (above,  below  and  on  both  sides,  with  wrapped 
edges)  (image  taken  from  [25]) 


4.1.3.  Comprehensive  learning  particle  swarm  optimization 

A  variant  of  the  PSO  called  the  comprehensive  learning  particle  swarm  optimization  uses  a 
novel  learning  strategy  in  which  all  other  particles’  historical  best  information  is  used  when 
updating  particle’s  velocity.  This  strategy  should  enable  the  diversity  of  the  swarm  to  be 
preserved  and,  in  the  end,  discourage  premature  convergence. 

There  are  numerous  proposed  variants  of  the  PSO,  but  almost  all  of  them  suffer  to  some 
extent  from  main  deficiency  of  the  PSO  -  premature  convergence  when  solving  multimodal 
problems.  The  reason  is  to  be  found  in  the  velocity  equation  for  the  classical  PSO.  All  the 
particles  learn  from  global  best  even  if  the  current  global  best  is  far  from  the  global  optimum. 
For  that  reasons,  the  particles  get  attracted  to  the  global  best  region  and  can  get  trapped  in 
local  optimum,  if  the  problem  to  be  optimized  is  complex  with  numerous  local  extremes.  One 
deficiency  of  the  classical  approach  is  outlined  in  [26]:  the  fitness  value  of  a  particle  is 
determined  by  all  N  parameters,  and  a  particle  that  has  discovered  the  region  corresponding  to 
the  global  optimum  in  some  dimensions  may  have  a  low  fitness  value  because  of  the  poor 
solutions  in  the  other  dimensions.  In  order  to  make  better  use  of  this  information,  the  authors 
have  proposed  a  new  learning  strategy  -  all  particles  personal  bests  are  used  to  update  the 
velocity  of  any  one  particle.  The  velocity  equation  is  changed  as  follows: 

vl  =  w-Vn+c-randO-{pbestfn‘(n) -x‘„)  (4.5) 

where  f,  =[fi(i),fi(2),...,fi(N)]  defines  which  particles’  personal  best  the  particle  i  should 
follow.  pbest{'(n)  can  be  the  corresponding  dimension  of  any  particle’s  personal  best  including 
its  own,  and  the  decision  depends  on  the  probability  Pc  called  the  learning  probability  which 
can  take  different  values  for  different  particles.  The  choice  is  made  as  follows.  For  each 
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dimension  of  each  particle  a  random  number  is  generated.  If  this  random  number  is  larger 
than  Pci,  the  corresponding  dimension  will  learn  from  its  own  personal  best,  otherwise  it  will 
learn  from  another  particle’s  best.  As  a  selection  procedure  for  the  case  the  random  number  is 
smaller,  tournament  selection  procedure  is  used.  Two  random  particles  (excluding  the  one 
whose  velocity  is  to  be  updated)  are  chosen  and  their  fitness  values  are  compared.  The  better 
particle  is  a  winner  and  its  personal  best  is  used  as  the  exemplar  to  leam  from  for  that 
dimension.  There  is  one  more  addition.  If  all  exemplars  of  a  particle  are  its  own  personal  best, 
one  dimension  is  randomly  chosen  to  leam  from  another  particle’s  personal  best.  In  order  to 
ensure  that  a  particle  learns  from  good  exemplars  and  to  minimize  the  time  wasted  on  poor 
directions,  we  allow  the  particle  to  leam  from  the  exemplars  until  the  particle  ceases  to 
improve  for  a  certain  number  of  generations  called  the  refreshing  gap  m  -  f  is  then  reassigned 
for  the  particle. 

According  to  authors  [26]  there  are  three  main  differences  between  the  CLPSO  and  the 
original  PSO.  Instead  of  using  only  personal  best  and  one  global  best  as  the  exemplars,  all 
particles’  personal  best  can  potentially  be  used.  Furthermore,  instead  of  learning  from  the 
same  exemplar  particle  for  all  dimensions,  each  dimension  of  a  particle  in  general  can  leam 
from  different  personal  bests  for  different  dimensions  for  a  few  generations.  Finally,  instead 
of  learning  from  two  exemplars  at  the  same  time,  each  dimension  of  a  particle  learns  from  just 
one  exemplar  for  a  few  generations.  These  three  changes  increase  the  swarm’s  diversity  and 
allow  improved  performance  when  solving  complex  multidimensional  problems.  This  is 
explained  as  follows.  In  classical  PSO,  if  the  personal  and  global  best  for  a  specific  dimension 
are  on  the  opposite  sides  of  the  particle’s  current  position  -  the  particle  may  oscillate. 
However,  the  global  best  is  more  likely  to  provide  a  larger  momentum  as  distance  between 
the  current  particle  position  and  global  best  position  is  likely  to  be  larger  than  distance  to 
personal  best.  From  these  reasoning  it  is  evident  that  global  best  may  influence  the  particle  to 
move  toward  it  even  if  it  is  in  a  local  optimum  region  that  can  be  far  from  the  global 
optimum.  If  both  personal  and  global  best  are  on  the  same  side  in  specific  dimension 
compared  to  particle  current  position,  the  particle  will  move  in  that  direction  and  it  may  be 
impossible  to  jump  out  of  the  local  optimum  area  once  particles  personal  best  falls  into  the 
same  local  optimum  region  where  the  global  best  is.  In  case  of  CLPSO,  the  particle  can  fly  in 
other  directions  by  learning  from  other  particles’  personal  best.  Furthermore,  as  each 
particle’s  personal  best  is  possibly  a  good  area,  the  CLPSO  is  neither  blind,  nor  random. 

CLPSO  introduces  one  new  parameter  called  learning  probability  Pc.  Question  naturally 
arises:  what  is  the  optimal  value  of  Pc.  Studies  have  showed  that  different  values  of  learning 
probability  yield  similar  results  on  simple  unimodal  problems  while  seriously  affecting  the 
CLPSO  performance  on  multimodal  problems.  In  order  to  address  the  problem  in  a  general 
manner,  learning  probability  is  typically  varied  between  particles.  In  this  way  particles  have 
different  levels  of  exploration  and  exploitation  ability  and  the  algorithm  is  more  suitable  for 
solving  diverse  problems.  The  following  expression  was  developed  in  order  to  set  the  Pci  for 
each  particle. 


Pa  =0.05  +  0.45  x 


f 10(z-l)^ 

\ 

-1 

exp 

l  ps- 1  J 

V 

J 

(exp(lO)-l) 


(4.6) 


In  this  way  the  learning  probability  between  particles  ranges  from  0.05  to  0.5  as  shown 
below. 
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Figure  4.5.  Learning  probability  Pc  for  each  particle  for  a  swarm  size  of  30. 


4.1.4.  Performance  comparison 

For  the  classical  PSO  we  have  used  parameters  very  close  to  one  recommended  in  [22], 
Boundaries  are  left  open  (invisible  walls),  and  no  fitness  is  calculated  outside.  The  inertial 
weight  w  decreases  linearly  from  0.9  to  0.4  during  one  optimization  run,  and  the  values  of  c/ 
and  C2  are  fixed  at  1.5.  In  contrast  to  the  recommendations  for  classical  PSO,  our  results 
suggest  that  increasing  the  particle  number  can  provide  better  results  since  the  search  space 
can  be  relatively  large.  Therefore,  in  the  optimizations  where  the  classical  PSO  has  been  used, 
the  population  size  was  typically  increased  to  at  least  100.  For  the  modified  classical  PSO 
with  local  best  topology  [25]  we  clamped  the  velocities.  The  local  neighborhood  was 
typically  set  at  5  and  was  defined  during  first  evaluation  based  on  normalized  distance 
between  particles.  In  the  case  of  CLPSO,  the  number  of  particles  was  keep  at  reasonable 
small  number  (around  40),  but  number  of  evaluations  was  set  at  much  higher  values  in  the 
range  of  30,000  or  more.  The  refreshing  gap  m  was  set  at  7  following  advice  given  in  [26], 
The  mentioned  steps  were  practical  to  implement,  as  the  developed  analyzing  kernels  were 
very  fast,  typically  needing  less  than  a  second  to  estimate  the  required  parameters  of  the 
structure.  It  should  be  noted  that  all  algorithms  have  been  able  to  obtain  similar  fitness  values 
for  most  cases,  so  it  is  reasonable  to  assume  that  local  extremes  have  largely  been  avoided. 
Details  about  typically  used  swarm  parameters  are  summarized  in  Table  4.1. 
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Classical  PSO 

Local  PSO 

CLPSO 

Parameters 

number  of 
particles 

o 

o 

II 

o 

o 

II 

o 

II 

maximum 

velocity 

V max*  none 

V  -05 

V max*  none 

minimum 

velocity 

V mini  none 

Vmin;  0.0001 

^ mini  none 

cognitive  ratio 

C]  =  1.5 

ci  =  1.5 

Ci  =  1.5 

social  ratio 

C2  =1.5 

c2  =1.5 

O 

II 

<N 

inertia  weight 

linearly 
decreasing  w: 

0.9  -►  0.4 

linearly 
decreasing  w: 

0.9  -►  0.4 

linearly 
decreasing  w: 

0.9  -►  0.4 

Properties 

neighborhood 

entire  swarm 

local  swarm  (5) 

entire  swarm 

solution  space 

[min,  max] 

ro.if 

[min,  max] 

boundary 

conditions 

invisible  wall 

invisible  wall 

invisible  wall 

Table  4.1.  Particle  swarm  optimization  parameters 


Fitness  (lower  is  better,  10  layer,  y  defines  losses) 


Classical  PSO 

Local  PSO 

CLPSO 

y  =  0 

0.009676 

0.009004 

0.009062 

y  =  0,00001 

0.416432 

0.314898 

0.314893 

y  =  0,001 

4.385527 

4.385527 

4.385527 

II 

o 

o 

9.409631 

8.527133 

9.409251 

Table  4.2.  Fitness  values  obtained  for  cloak  with  10-layers  with  various  losses 


Fitness  (lower  is  better) 

No.  of  layers 

Local  PSO 

CLPSO 

7 

1.163527 

1.168191 

8 

0.551858 

0.551849 

9 

0.328954 

0.328936 

10 

0.009004 

0.009062 

12 

0.000966 

0.001020 

Table  4.3.  Fitness  values  obtained  for  lossless  cloak  with  different  number  of  layers 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Metamaterial-Based  Cylinders  Used  for  Invisible  Cloak  Realization 


41 


4.2  Optimized  cloak 

Recently,  several  research  groups  investigated  optimization  of  cylindrical  cloaks. 
However,  they  have  used  either  fully  anisotropic  materials  [27]  (with  no  indication  how  to 
practically  realize  it  in  microwave  frequency  range),  or  they  have  optimized  cloaks  for  optical 
or  infrared  range  of  frequencies  [28]  where  the  permittivity  smaller  than  one  can  be  obtained 
using  plasmonic  materials  (e.g.  isotropic  materials  can  be  obtained  using  the  noble  metals, 
polar  dielectrics  or  some  semiconductors)  or  composite  materials  obtained  by  properly 
embedding  metallic  nanoparticles  and  nanowires  into  dielectric  background.  However,  in 
microwave  frequency  region  the  obtained  structures  with  permittivity  smaller  than  one  are 
anisotropic.  Therefore,  we  have  concentrated  our  investigation  on  the  cloaks  in  microwave 
frequency  region  that  can  be  obtained  relatively  easily  using  periodic  structures  that  change 
constitutive  parameters  in  only  one  direction.  All  the  presented  results  are  for  Schurig-type  of 
cloaks,  i.e.  for  cloaks  that  have  only  radial  component  of  permeability  that  varies  along  the 
radial  direction.  Similar  results  can  be  obtained  for  TE  cloaks  with  only  radial  component  of 
permittivity  that  is  a  function  of  radial  coordinate. 

Using  the  developed  PSO  algorithm  we  have  optimized  the  relevant  constitutive 
parameters  (jur  and  ez)  of  each  layer  of  the  circular-cylindrical  cloak,  in  order  to  obtain  a  cloak 
that  does  not  scatter  the  electromagnetic  field  at  the  central  frequency,  i.e.  the  “perfect”  cloak 
at  normal  incidence.  It  should  be  noted  that  such  a  “perfect”  cloak  is  designed  with  the 
assumption  of  lossless  metamaterial  {y  =  0).  The  calculated  relevant  constitutive  parameters 
of  the  Schurig  cloak  and  the  optimized  cloak  are  given  in  Table  4.4.  The  considered  number 
of  layers  is  10,  i.e.  both  the  Schurig  cloak  and  the  optimized  cloak  consist  of  the  same  number 
of  layers.  The  layer  no.  1  represents  the  innermost  layer  while  the  layer  no.  10  is  the 
outermost  one.  The  level  of  complexity  of  the  construction  of  the  optimized  cloak  is  the  same 
as  in  the  construction  of  the  Schurig  cloak  (if  needed,  the  interval  of  obtainable  permeability 
values  can  be  easily  adjusted  in  the  optimization  process  if  some  of  the  values  of  jur  cannot  be 
obtained  in  the  practical  cloak  realization).  However,  the  physical  insight  into  wave 
propagation  through  the  optimized  cloak  is  lost.  The  calculated  total  scattering  width  of  the 
optimized  cloak,  together  with  the  comparison  to  the  original  Schurig  cloak,  is  given  in  Fig. 
4.6.  The  results  are  given  in  a  narrow  frequency  range  above  the  central  frequency  (8.5  GHz), 
since  the  cloak  behaves  quite  symmetrically  for  frequencies  lower  than  the  central  one.  The 
invisibility  gain  is  simply  obtained  by  dividing  the  total  scattering  widths  of  the  referent  case 
(PEC  only)  and  the  case  with  the  cloak  present.  It  is  evident  that  by  using  the  optimization 
process  it  is  possible  to  design  a  virtually  invisible  cloak  with  the  obtained  invisibility  gain  of 
about  1,400  (although  the  physical  insight  of  the  cloak  is  lost).  Note  that  the  bandwidth  of  the 
optimized  cloak  is  reduced  even  more,  compared  to  the  Schurig  cloak.  While  in  the  Schurig 
cloak  the  estimated  invisibility  bandwidth  is  approximately  0.2%,  for  the  optimized  cloak  the 
bandwidth  is  only  about  0.08%. 

The  angular  variation  of  the  bistatic  scattering  width  ( 02D )  has  been  calculated  at  the 
central  frequency  just  to  check  if  there  is  any  direction  with  distinctively  stronger  scattered 
field  (Fig.  4.7).  As  can  be  seen,  the  scattered  field  of  the  optimized  cloak  is  around  20  dB 
lower  than  that  of  the  Schurig  cloak. 
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Layer 

No. 

Schurig  cloak 

Optimized  cloak 

Optimized  cloak 
(y  =  0.0003) 

=3.431; 

£z  =3.858; 

fz  =3.813; 

A*  =1-0 

A>=l-0 

Mr 

Mr 

Mr 

1 

0.003 

0.056 

0.020 

2 

0.022 

0.025 

0.016 

3 

0.051 

0.010 

0.020 

4 

0.085 

0.092 

0.078 

5 

0.119 

0.587 

0.192 

6 

0.154 

0.046 

0.068 

7 

0.187 

0.216 

0.193 

8 

0.219 

0.999 

0.500 

9 

0.249 

0.076 

0.093 

10 

0.278 

0.616 

0.503 

Table  4.4.  Relevant  constitutive  parameters  of  the  Schurig  cloak  and  the  optimized  cloaks 


Frequency  (GHz) 

Figure  4.6.  Normalized  total  scattering  width  vs.  frequency  for  the 
Schurig  cloak  and  the  invisibility  gain  -  optimized  cloak.  The  constitutive 
parameters  of  both  cloaks  are  given  in  Table  4.4 
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Figure  4.7.  Normalized  bistatic  scattering  width  of  the  considered  cloak 
realizations  at  the  central  frequency  (8.5  GHz).  The  constitutive 
parameters  of  both  cloaks  are  given  in  Table  4.4 


4.2.1.  Number  of  layers 

The  Schurig  cloak  consists  of  ten  metamaterial  layers,  which  is  the  starting  point  of  our 
study.  Here  we  shall  study  how  the  number  of  layers  in  the  cloak,  i.e.  the  subtlety  of  the 
approximation  of  radial  permeability,  relates  to  the  level  of  invisibility  obtained  by 
optimization.  All  the  results  have  been  obtained  by  PSO  optimization,  keeping  the  total 
dimensions  of  the  cloak  intact  (i.e.  the  inner  and  outer  radii).  From  Fig.  4.8  it  can  be  seen  that 
the  normalized  total  scattering  width  decreases  as  the  number  of  layers  increases.  It  is 
interesting  to  note  that  better  invisibility  performance  than  with  the  Schurig  cloak  (that  has 
<y T  j2ci  =  0.717)  can  be  obtained  with  only  three  metamaterial  layers.  Also,  it  should  be 
stressed  that  the  level  of  achieved  invisibility  is  significantly  increased  once  the  tenth  layer  is 
added.  With  ten  layers  the  cloak  becomes  practically  invisible,  although  only  at  a  single 
frequency. 
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Number  of  layers 

Figure  4.8.  Normalized  total  scattering  width  vs.  number  of  layers  for  a 
cloaked  cylinder  with  1.5  X  diameter.  The  parameters  of  layers  where 
optimized  with  PSO  at  a  single  frequency  of  8.5  GHz. 


We  have  also  investigated  how  the  radius  of  the  cylinder  that  we  would  like  to  cloak 
influences  the  needed  number  of  cloak  layers.  The  considered  cylinder  in  Schurig  case  had  ~ 
1.5  X  diameter.  In  Fig.  4.9  we  have  investigated  the  cylinder  with  1.0  X  diameter  (at  central 
frequency  8.5  GHz),  while  the  ratio  of  the  outer  and  inner  cloak  radius  was  kept  the  same  (i.e. 
the  considered  radii  are  a  =  1.765  cm  and  b  =  3.836  cm).  It  is  interesting  to  note  that  with  the 
reduction  of  the  cloaked  cylinder  radius  the  needed  number  of  layers  in  order  to  obtain  “the 
perfect  cloak”  is  also  reduced,  and  in  this  case  one  needs  only  7  layers  to  obtain  invisibility 
gain  more  than  1000. 
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Number  of  layers 

Figure  4.9.  Normalized  total  scattering  width  vs.  number  of  layers  for  a 
cloaked  cylinder  with  1.0  X  diameter.  The  parameters  of  layers  where 
optimized  with  PSO  at  a  single  frequency  of  8.5  GHz. 


4.2.2.  Losses  and  tolerance 

In  addition,  we  have  estimated  the  influence  of  losses  in  the  material  on  the  maximum 
achievable  invisibility  gain.  We  have  started  from  the  lossless  cloak  and  gradually  increased 
losses  running  the  optimization  procedure  for  all  the  considered  values  of  y  independently. 
The  losses  are  described  with  factor  yin 


^eff  “I 


f2-f02-jrf' 


(4.7) 


As  expected,  the  total  scattering  width  increases  (and  therefore  the  obtained  invisibility  gain 
decreases)  as  the  losses  increase  (see  Fig.  4.10).  It  has  also  been  found  that  the  losses  of 
Y=  io  3/0  are  acceptable,  for  the  invisibility  gain  of  the  optimized  lossy  cloak  has  fallen  to 

approximately  the  same  value  as  the  invisibility  gain  of  the  lossless  Schurig  cloak  (around  3). 
In  the  same  figure  we  have  plotted  the  influence  of  material  losses  on  the  bandwidth  of  the 
cloak  (i.e.  the  frequency  range  where  invisibility  gain  is  larger  than  1).  It  can  be  seen  that  the 
bandwidth  increases  with  the  introduction  of  the  losses,  while  at  the  same  time  the  maximum 
invisibility  gain  is  reduced. 


Furthermore,  sensitivity  of  the  optimal  solution  to  construction  tolerances  has  been  tested. 
We  have  run  statistical  analysis  to  estimate  how  much  the  variation  (as  a  result  of 
construction  tolerances)  of  the  values  of  the  radial  component  of  permeability  tensor  /u,  and  of 
permittivity  sz  influences  the  obtained  total  scattering  width.  Fig.  4.11  shows  the  mean  value 
of  total  scattering  width  if  /ir  and  ez  are  varied  randomly  by  0%  to  5%  of  their  initial  (ideal) 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Metamaterial-Based  Cylinders  Used  for  Invisible  Cloak  Realization 


46 


value  with  uniform  distribution  (worst  case).  Two  cases  are  shown  -  with  and  without  losses. 
It  can  be  seen  that  the  cloak  is  not  too  sensitive  to  parameter  tolerances  and  lossy  cloaks  are 
even  less  sensitive.  We  can  conclude  that  cloaks  that  are  quite  invisible  can  be  realized  even 
with  reasonably  large  parameter  variances.  The  cloak  with  small  losses  is  also  easier  to 
realize  in  practice  since  the  variation  of  parameters  is  much  smaller,  i.e.  all  the  values  of 
radial  permeability  lie  in  range  from  0  to  0.5  (Table  4.4,  column  3),  so  the  design  could 
proceed  as  in  [14]. 
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Figure  4.10.  Dependence  of  obtained  total  scattering  width  and  bandwidth 
on  losses  for  the  invisibility  gain  -  optimized  cloak. 


Figure  4.11.  Dependence  of  obtained  total  scattering  width  on 
construction  tolerances  for  the  invisibility  gain  -  optimized  cloak.  The 
mean  value  of  total  scattering  width  is  shown. 
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4.2.3.  Bandwidth  -  optimized  cloak 

We  have  also  investigated  the  possibility  of  enlarging  the  bandwidth  of  the  cloak  while 
keeping  the  same  level  of  the  invisibility  gain  (i.e.  the  minimum  of  the  total  scattering  width) 
as  the  one  of  the  Schurig  cloak.  Two  cases  have  been  considered  -  the  optimization  of  the 
lossless  cloak  (/=  0)  and  the  optimization  of  the  lossy  cloak  with  small  losses  of  y=  10"3/o  , 

which  have  been  shown  to  be  still  acceptable.  The  calculated  total  scattering  widths  are 
shown  in  Figs.  4.12  and  4.13,  respectively.  The  comparison  with  the  Schurig  cloak  (lossless 
and  lossy,  respectively)  is  also  given.  The  optimized  bandwidth  is  around  0.4%  for  the 
lossless  case  and  about  0.58%  for  the  lossy  case.  Meanwhile,  the  Schurig  cloak  exhibits 
bandwidths  of  0.2%  (lossless  case)  and  0.24%  (lossy  case).  In  other  words,  the  improvement 
in  the  bandwidth  compared  to  the  Schurig  cloak  is  between  2  and  2.5  in  both  considered 
cases. 


Figure  4.12.  Normalized  total  scattering  width  vs.  frequency  for  the  Schurig 
cloak  and  the  bandwidth  optimized  cloak.  Lossless  cloaks  are  considered  ( y  =  0). 
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Figure  4.13.  Normalized  total  scattering  width  vs.  frequency  for  the  Schurig 
cloak  and  the  bandwidth  optimized  cloak.  Lossy  cloaks  are  considered 

(r=iO'3/0  )• 


4.2.4.  Three  layers  cloak 

Fig.  4.8  shows  the  minimum  total  scattering  width  that  is  possible  to  obtain  by  respective 
number  of  metamaterial  layers  (the  goal  of  optimization  in  that  optimization).  It  can  be  seen 
that  the  normalized  total  scattering  width  decreases  as  the  number  of  layers  increases.  From 
section  4.2.1  it  also  follows  that  by  using  optimization  it  is  possible  to  achieve  better 
invisibility  performance  than  with  the  Schurig  cloak  (that  has  aT  /2 a=  0.717)  by  employing 
only  three  metamaterial  layers.  As  3-layers  cloak  is  undoubtedly  more  convenient  to  produce 
in  practice,  it  is  obviously  interesting  to  provide  study  of  the  three  layer  cloak  design  and 
performance  in  more  details.  The  results  of  the  three  layer  optimization  are  summarized  in 
Table  4.5. 


Layer 

No. 

Optimized 

cloak 

=3.841; 

A,  =1-0 

Mr 

1 

0.021 

2 

0.123 

3 

0.218 

Table  4.5.  Relevant  constitutive  parameters  of  the  3 -layer  cloak 
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The  obtained  relative  permeability  values  are  reasonably  easy  to  obtain  in  practice,  making 
realization  of  such  cloak  possible  by  means  of  split  ring  resonators.  The  dispersion  of  3-layers 
cloak  was  also  taken  into  account  and  modeled  by  Lorentz  model.  In  Fig.  4.14  frequency 
dependence  of  total  scattering  width  of  3 -layers  cloak  is  given,  as  well  as  the  comparison  with 
Schurig  cloak  (which  is  ten-layered).  The  results  are  given  in  a  narrow  frequency  range  above 
the  central  frequency  (8.5  GHz),  since  the  cloak  behaves  quite  symmetrically  for  frequencies 
lower  than  the  central  one.  Obviously  by  significantly  smaller  number  of  layers  similar  results 
have  been  obtained  (the  invisibility  gain  was  slightly  improved  at  the  expense  of  bandwidth). 


Figure  4.14.  Normalized  total  scattering  width  vs.  frequency  for  the 
Schurig  cloak  and  the  invisibility  gain  -  optimized  3-layer  cloak 


Furthermore,  sensitivity  of  the  optimal  solution  to  construction  tolerances  has  been  tested.  We 
have  run  statistical  analysis  to  estimate  how  much  the  variation  (as  a  result  of  construction 
tolerances)  of  the  values  of  the  radial  component  of  permeability  tensor  jur  and  of  permittivity 
ez  influences  the  obtained  total  scattering  width.  Fig.  4.15  shows  the  mean  value  of  total 
scattering  width  if  ju,  and  sz  are  varied  randomly  from  0%  to  5%  of  their  initial  (ideal)  value 
with  uniform  distribution  (worst  case). 
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Figure  4.15.  Dependence  of  obtained  total  scattering  width  on  construction 
tolerances  for  the  3-layers  cloak.  The  mean  value  of  total  scattering  width  is  shown. 

As  an  example,  2%  tolerance  was  studied  in  more  details.  Constitutive  parameters  were 
varied  uniformly  in  2%  interval  of  exact  values  for  106  times  and  total  scattering  width  at 
central  frequency  (8.5  GHz)  was  calculated  for  each  case.  The  number  of  values  per  small 
interval  of  total  scattering  width  is  given  in  Fig.  4.16.  Obviously  most  values  are  concentrated 
around  the  exact  value  (around  0.47),  suggesting  that  with  reasonable  construction  tolerances 
it  is  still  possible  to  obtain  good  invisibility  performance. 


Figure  4.16.  Distribution  of  106  total  scattering  widths  obtained  by  2% 
variation  of  constitutive  parameters 
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We  have  also  investigated  the  influence  of  the  cloaked  cylinder  diameter  on  performance 
of  the  cloak  with  three-layers  (we  have  kept  the  ratio  b/a  =  2.173  as  in  the  Schurig  cloak).  The 
results  are  summarized  in  table  4.6.  It  can  be  seen  that  the  quality  of  the  cloak  with  three- 
layers  is  drastically  improved  if  we  reduce  the  diameter  of  the  cloaked  cylinder:  for  the 
diameter  of  the  cloaked  cylinder  2 a  =  0.5  X  the  total  scattered  width  is  crT  /  2 a  =  0.053 ,  while 
for  the  2 a  =  3.0  X  the  total  scattered  width  is  aT /2 a  =  0.876  (this  variation  is  larger  than  10 
times).  On  the  other  side  the  total  scattered  width  of  the  10-layer  Schurig  cloak  varies  only 
between  0.939  and  0.624.  Furthermore,  we  have  investigated  how  many  layers  the  optimized 
cloak  needs  to  contain  in  order  to  have  better  performance  than  the  10-layer  Schurig  cloak. 
We  have  concluded  that  for  a  small  diameter  of  the  cloaked  cylinder  it  is  enough  to  have  2 
layers,  while  for  larger  diameters  more  layers  are  needed  (e.g.  for  3  X  cloaked  cylinder  one 
need  to  have  6-layer  structure). 


Cloaked  cylinder 
diameter  (2  a) 

oT  /2  a 

of  optimized  cloak 
with  3  layers 

<jt  jla 

of  Schurig  cloak 
with  10  layers 

Minimum  number 
of  layers 

of  optimized  cloak 

0.5  X 

0.053 

0.939 

2 

1.0  X 

0.397 

0.716 

3 

1.5  X 

0.469 

0.717 

3 

2.0  X 

0.745 

0.703 

4 

3.0  X 

0.876 

0.624 

6 

Table  4.5.  Total  scattering  widths  of  optimized  3-layer  cloaks  as  a  function  of  cloaked 

cylinder  diameter. 

Finally,  we  have  investigated  what  is  the  optimal  size  (outer  radius)  of  the  cloak.  As  an 
example  we  have  considered  a  PEC  cylinder  of  variable  radius  and  a  three-layer  cloak.  The 
variables  to  be  optimized  are  permeability  in  radial  direction  (gq,  =  1.0),  the  permittivity  sz  and 
the  ratio  b/a.  The  results  are  given  in  table  4.6.  It  can  be  seen  that  for  small  radius  of  the 
cloaked  cylinder  the  cloak  has  resonant  nature  (i.e.  the  values  of  permeability  in  radial 
direction  are  alternating)  and  the  optimal  cloaks  have  large  permittivity  values.  On  the  other 
side,  for  larger  sizes  of  the  cloaked  cylinder  the  values  of  permeability  in  radial  direction 
smoothly  grows  with  the  radial  coordinate,  and  the  outer  radius  and  permittivity  values  are 
quite  similar  to  the  Schurig  design  ( sz  =3.43 1 ,  b/a  =  2. 1 7). 


Cloaked 

cylinder 

diameter 

(2a) 

(7  j*  j2  a 

of 

optimized 
cloak  with 

3  layers 

b/a 

outer 

radius 

hr 

of  the 
first 
layer 

hr 

of  the 
second 
layer 

hr 

of  the 
third 
layer 

<jt  /2  a 

of  Schurig 
cloak  with 
10  layers 

0.5  X 

0.009 

1.924 

0.8687 

0.0050 

0.3186 

8.4627 

0.939 

1.0  X 

0.206 

1.580 

1.0 

0.0145 

0.1566 

8.0124 

0.716 

1.5  X 

0.440 

2.246 

0.0251 

0.1332 

0.2449 

3.5321 

0.717 

2.0  X 

0.561 

1.848 

0.0242 

0.1071 

0.1822 

4.5858 

0.703 

Table  4.6.  Relevant  constitutive  parameters  of  the  3 -layer  cloak  with  variable  outer  radius 
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5  OBLIQUE  INCIDENCE  PLANE  WAVE  SCATTERING  BY 

SCHURIG  CLOAK 
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Oblique  incidence  plane  wave  scattering  by  Schurig  cloak 


5.1.  ANALYSIS  METHOD 

In  the  preceding  chapters  mathematical  description,  analysis  and  optimization  of  a 
cylindrical  cloak  for  perpendicular  incidence  is  made.  In  this  chapter  we  will  analyze  the  case 
when  the  incoming  plane  wave  has  an  arbitrary  angle  of  arrival  (Figure  5.1).  Oblique 
incidence  is  more  complicated  to  analyze  than  perpendicular  incidence  because  TEZ  or  TMZ 
incoming  wave  will  also  induce  cross-polarization  component.  Because  of  this  the  resulting 
differential  equations  will  be  coupled.  We  would  like  to  comment  that  this  is  a  still  open 
problem  in  scientific  literature  -  till  now  only  J.  Monzon  tried  to  solve  this  complex 
electromagnetic  problem.  In  [29]  he  derived  the  solution  using  Bames-type  representation  of 
cylinder  functions  (with  the  main  practical  question  opened  -  how  to  make  an  efficient 
program  that  will  solve  this  complex  electromagnetic  problem). 


Figure  5.1  A  sketch  of  the  analyzed  structure 


Our  analysis  begins  with  the  Maxwell  equations.  For  the  analysis  only  the  Faraday  and 
Ampere  law  are  used: 

V  x  E  =  -  j  (5.1) 

V  x  H  =  j  ry  £  E  (5.2) 

In  equations  (5.1. a)  and  (5.1.b)  and  £  represents  tensor  of  permeability  and  permittivity  with 
3  components. 
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PPP 

o 

0 


0 
0 

o  /4z 


£  = 


PP 

0 


0  0 

SM>  0 
0  s„ 


(5.3) 


(5.4) 


Cloak  has  cylindrical  symmetry  and  because  of  that  simplest  way  to  solve  the  problem  is  to 
use  the  cylindrical  coordinate  system.  Using  the  definition  of  the  V  x  operator  in  the 
cylindrical  coordinate  system  [30]  we  get  from  equation  (5.1): 


1  dEz  dE<f>  . 

p  St  -  & 

(5.5) 

8EP  aE, 

&  -  ep  -J Wrt 

(5.6) 

1  8(pEf)  1  dE 
p  op  p  C(j) 

(5.7) 

and  from  equation  (5.2): 

1  5Hz  dH . 

P  <¥  ’  az  = 

(5.8) 

dHp  SH ,  _ 

&  ap 

(5.9) 

1  «K) 

~  A  A  ~  J^O^ZZ^Z 

p  op  p  C(j) 

(5.10) 

From  the  previous  equations  we  can  express  the  p  and  (j)  of  the  E  and  H 

fields  in  the 

following  way: 

r  .  1  1  dHz  (  ,  1  dH 

Ep  =  -J - -  +  j - 5:1 

P  d<P  C0£App  dz 

(5.11) 

r  •  1  dHp  ,  ■  1  dHz 

dz  c O£0£ H  dp 

(5.12) 

u  .  1  1  dEz  .  1  dE 

°>PoPpP  P  d<f>  ojpX)p  dz 

(5.13) 

Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Metamaterial-Based  Cylinders  Used  for  Invisible  Cloak  Realization 


55 


H4=  j 


1  dEn 


-J- 


1  dE 


dz  ®PoPw  8P 


(5.14) 


Figure  5.2.  Translation  symmetry  in  z-direction 


From  the  translation  symmetry  in  z-direction  we  can  assume  dependence  of  E  and  H  fields  in 
z  coordinate  as: 


(5.15) 

H(pJ,z)  =  TL(pJ)e-^ 

(5.16) 

Symbol  E  and  H  represents  E  and  H  fields  without  dependency  on  z  coordinate.  With  that 
assumption  we  get  equations: 


E  -  j  '  18H-  ,  k>  H 

P  <*>e0 spp  p  d(j)  co£{)£pp  * 

(5.17) 

p  _  ~kZ  TT  |  •  1  dH, 

CO£{)£h>  dp 

(5.18) 

H  =  j  1  1  d^z  k=  E 

m>pPP  p  8</>  m>pPP 

(5.19) 

P  kz  d  ■  1 

Ph  ®PoPm  8P 

(5.20) 

Because  we  want  to  represent  the  p  and  (j)  components  of  E  and  H  fields  as  a  function  of 
E_  and  H_ ,  we  will  substitute  equations  (5.17)  and  (5.18)  into  the  equations  (5.19)  and  (5.20) 
and  vice  versa.  After  short  manipulation  we  get: 
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£  ■  °EkPH,  1  8Hz  kz  8Ez 

P  ml-k*  p  8<j>  m)p-k]dp 

(5.21) 

p  _  :  k-  1  8Ez  cop0ppp  8Hz 

*  n£P<rk]  P  d<j>  m^-k2  dp 

(5.22) 

tt  ■  m£*£H.  1  dEz  kz  8Hz 

p  P  a*  m2p,-k2  8p 

(5.23) 

TT  _  ■  K  1  8Hz  (0S,epp  8Ez 

*  m]p-k;  P  8(p  m)p-k2z  dp 

(5.24) 

where  the  used  abbreviations  are: 

(5.25) 

m]P  =  (o1PpPH’£»£pp 

(5.26) 

Equations  (5.21)  -  (5.24)  are  substituted  into  equations  (5.11)  and  (5.14)  and  after  short 
derivation  we  get  two  coupled  partial  differential  equations 


P 


P 


0  /v  /y 

82E  8E 

+  p-J-  + 


dp 

2  /  2 

-k 


dp 

k 


p2(™lP-k;) 


E+ — $ 


-k]  su  82E, 


pp 


8 


mpS~K  C0£iPpp  $</> 


8H z 

dp 


+  P 


8 


■K  spp  8(j)2 


G>£0£pp  dp 


8EL 

8(f) 


=  0 


(5.27) 


2  82Hz  8Hz  2  /  2  1  2  \  Pzz  tt  ~k 

p^r+p^z+p(ml>-k--)£rH>  + 


8  H, 


m 

p  2 


dp2 
2  k: 


dp 

k_ 


M. 


pp 


m. 


K  ppp  d</> 


-  + 


8 


Wfip-K  cop0ppp  8<t> 


8EZ 

dp 


P 


8 


WPPpp  dp 


8E_ 
Kdt  j 


=  0 


(5.28) 


Using  the  Schwartz  theorem  from  multivariable  calculus  [31]  and  assumption  that  solutions  of 
above  equations  are  from  C2  ((a,oo)x(0,27i))  ,  where  a  represents  the  inner  radius  of  cloak, 

we  can  interchange  the  order  of  derivatives  and  get: 


P 


2  82Ez  8E  2/2  ,  2  \ 

w+p^+p 


s  -  m, 
E  +-  * 


-  K 2  £u  82E, 


PP 


mp*~K  £PP  d<t> 


P- 


mlP~k2z  kz  8 

{  dHz  1 

k  8 

f  8Hz "1 

™2pt  ~k:  cos ,spp  dp 

l^J 

1  P  Z 

cos0spp  dp 

(5.29) 


=  0 
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,2  5  Hz  d. Hz  2  /  2  7  2  \  P-_-_  rV  . 


+^+^K-i-)^=+ 


mlp~kl  MPO  df 


K*-kl 

K 

d  ( 

mlp~kl 

VPoPpp 

dp  v 

d  dE, 


(5.30) 


Now  we  have  2  coupled  partial  differential  equations  which  are  very  hard  to  solve.  We  want 
somehow  to  remove  derivation  dependent  on  <j)  and  get  2  coupled  ordinary  differential 
equations.  From  the  spectral  representation  method  of  linear  operators  [32],  [33]  and 
theorems  from  functional  analysis  [34]  we  know  that  function  can  be  decomposed  in  the  base 
of  eigenvectors  of  the  operator. 


After  we  substitute  (5.31)  and  (5.32)  into  eqr 

p\2  f  +oo  \  rs  f  +°° 

p2-^t  Y k{pyn^  +p—  yk (p) 
do2  i  At  ~z  ’  do\  At  ~z  y  } 


2j 

dp  y„=-x 


+00  \ 

Xtipp*  \+p 


a  ( s 


k(p,t)=1Lk(p)j* 

n=-  oo 

+00 

n=—  oo 

52)  into  equations  (5.29)  and  (5.30)  we  get: 


(5.31) 

(5.32) 


2  /  2 

P  K 


;  mpt~E  £pp 

%(p) 

~K  <*>£o£PP  dp\  d(!)\n=_Vj 


m2,n  -  k2  s , ,  d2  (  4^  -  /  x 

Is  A 

m  ,-k  s  do  l  „tt. 


(5.33) 


P  YH _ ^ _ Z 

Spp 


k  r>  (  d  f  42.  - 

+/0_s — —  —  y //;(/?)< 


,j"^  _  Q 


Pp  k  +00  \  P  f  +GO 

p2ir  +PT-  ZC(p) 

V«=-oo  y  v«=-°° 

f  H-A o)e'A+Eetd±, 


//  +Q0 
M' pp  \n=-cc 

Kt~kl  K  d(df^ 

p— y - j - 2 - —  > 

m,i>P  ~ k-_  ®A0 PPP  dp\d<t>\n=- 


(5.34) 


->)<?*  \+nY-Y—^T\  Y  Hnz{p)e^  +  (5.34) 

)  ™<t)p  ^ z  Ppp  dfi  \n=-oo  J 

+°°  l  p  f  p  f  +qo 

ls(/+w  -p— s-f-  t-  zt+)7w  =o 

n=-x>  JJ  ^PoPpp  dp  V+ ^  V  7?=-CO  /  y 

^dependent  and  operators  with  derivatives  are  linear  [35]  so  we 
rivatives,  and  by  using  linear  independency  we  get: 


Eigenfunctions  are  lineary  independent  and  operators  with  derivatives  are  lint 
can  interchange  sums  and  derivatives,  and  by  using  linear  independency  we  get: 

2d2t  .  dt  .  (  u  2  l2^zz  2m]p-kl  £^)dn  . 

-^-n2 — ^ + 

_  „  2  7  2  _  ~z 


2  d2En  dEn 

P  -rr+P-f- 
dp  dp 


2K~kl) 


~kl  £n 


}"k-_  Yk 


p\  — - 2 - ^ - - - - - l 

\°>e*epp  m^~k 2  (O£0£, 


~k-  k  Xw:  0 


(5.35) 
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2  d2H"  dH " 

p — =^  +  p^^-+ 
dp  dp 

■2 


mi-k: 


P2(m^-k')  —  -n  2  ,2 

Ppp  m^~K  P, 


p 


f  Kt-K  jnkz 


]nkz 


mln~k:  m>Pnn  ®P0P, 


dE" 


pp  J 


dp 


=  0 


m  + 


z  e~pp  j 


In  order  to  get  compact  notations  we  introduce  the  following  coefficients: 

,  /  \  2  /  2  i  2  \  ^zz  -  _  k- 

ae(p)  =  p \%p-K)— 


pp 


m*~k$  £pp 


Ah  (p)  =  P 


f  2  2  ^ 

Mz  m^-kz  Mz 

\  C0£<Ppp  mp0  —  K  coeospP  j 


Poo  '%>  ~K  P, 


Z  r~pp 


Be  (p)  =  P 


]nk- 

mlP  ~  kl  tvPoPpp  ®PoP, 


j  nk2 


~pp  j 


Where  line  above  coefficient  represents  function  in  the  continuous  domain.  The 
notation  of  the  considered  partial  differential  equations  is: 


P 


2  d2En  dE:  -  a-  -  dH: 


+  P^r~  +  aeEhz  +  ah 


dp  dp 


dp 


=  0 


d2Hn 


P 


-  +  P- 


dHn,  - 


+bhh:+be 


=  dE: 


dp2  'dp  n  dp 

Outside  cloak  ji  =  1  and  £  =  1  which  simplifies  equations  to: 


=  0 


where: 


p 


dp  dp 


d2Hn  dH" 

+  p^  + 


dp  dp 


K  = 


2k 

T 


For  TMZ  incident  wave  we  can  write  (see  Figure  5.3) 

Ef  =  E'(j  (xcos(a)  +  2sin(a))eAkoSin(a)x+ik°cos(a)z 


(5.36) 

(5.37) 

(5.38) 

(5.39) 

(5.40) 
compact 

(5.41) 

(5.42) 

(5.43) 

(5.44) 

(5.49) 

(5.50) 
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The  scattered  field  is  equal 


(5.51) 


Eo=R-K 

E‘v  =  E;  (-icos (a)  +  zsm(a))e^sinia)x+ik°cosia)z 


(5.52) 


Figure  5.3.  Oblique  incidence  on  cloak,  TMZ  mode 


We  are  interested  in  E_  and  H_  components. 


E[  =  El0  sin(a)rf  j*°sin(“)x+j*oCOs(“)z 

(5.53) 

El  =  El  sin(a)ej*°sin(“)x+j*°cos {a)z 

(5.54) 

We  can  choose  z=0  for  incidence  plane 

Ei=El0sm(a)e-ik°sm(a)x 

(5.55) 

El  =  El  sin (cc)e]k^Ha)x 

(5.56) 

Solutions  for  equations  (5.43)  and  (5.43)  are  [36] 

&  =  J,(Kp) 

(5.57) 

e:=r,»?(Kp) 

(5.58) 
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where  (5.57)  represents  incoming  and  (5.58)  outgoing  wave.  From  representation  exponential 
functions  with  Bessel  function  we  get  [37]: 

oo 

E;  =  E„  sin(o)  X  J-  +  (5.59.a) 

n=— oo 

00 

Hnz  =  E0  sin(a)  £  (5.59.b) 


Figure  5.4.  Oblique  incidence  on  cloak,  TEZ  mode 


For  TEZ  polarization  similar  procedure  gives  us: 

oo 

K  =  E0  sin  (a)  X  j'XH ?(kzp)e]n*  (5.60.a) 

n=-  oo 

00 

Hnz  =E0sm(a)Xr(j„(KP)  +  Cnn(n2)(kzp))ein*  (5.60.b) 


Equations  inside  the  cloak  region  will  be  solved  using  finite  difference  method  [38],  [39], 
Derivation  is  approximated  with  central  difference  formula 


dp  ^°  ~  2 h 


(5.61) 


A 2|;  [a ] + 1;  [ a  - *] 


lr 


(5.62) 


where  h  represent  difference  between  two  points  in  the  radial  direction  in  equidistant  domain 
subdivision.  Distance  from  PEC  is  approximated  with: 
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p  =  a  +  ih,  i  =  0, 1,  2, . . . 

Coefficients  in  equations  (5.37)  -  (5.40)  are  approximated  with: 


Ae  [/]  =  (a  +  ihf  (ml  ~k2)^~ 

bpp  1p4'  z  bpp 


Ah  [z]  =  (a  +  ih) 


)nkz  m)p-k]  j  nkz 


\CO€q£pp  mp(p  k_  CQS^Spp  j 


B„  \i]  =  {a  +  ih)2(m2p<p-kz) 


i\  /4 


-n 


m2,  -  k2 

2  P&  z 


M, 


PP 


m*P-K  bpp 
\ 


f  ml~k2  j  nkz  ]nk. 


BE[i]  =  (a  +  ih) 

\mtp-K  °BkMPP  ojp()ppp  J 

By  this  equations  (5.41)  and  (5.42)  are  transformed  into: 


h 1 


2  h 


q,]|:[i]tq,]lfcrtN=o 


2  h 


2  h 


r  n  ,  r  n  r  n£"f/  +  ll-£"f/-ll 
bh  [/] h;  [/] + be  [/]^i — J..-‘L  J  =  o 


2  h 


(5.63) 

(5.64) 

(5.65) 

(5.66) 

(5.67) 

(5.68) 

(5.69) 


By  grouping  the  coefficients  we  get: 


KM- 

£  [■•+!]• 


a  +  ih 


V 

r  a  +  ih 


a  +  ih 


h  , 

\2 


h 


+  - 


2  h  j 

a  +  ih 
~2h~ 


}  +  E:[i]\-2 


a  +  ih 


+  Ae  [/]  [  + 


(5.70) 


2  h 


2  h 


=  0 
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CM- 


r  a  +  ih 
h  j 

a  +  ih 

v  h  j 


+  - 


a  +  ih 
~2. h 

a  +  ih 
~2h 


+  Hnz\i]\-2 


a  +  ih 


\  n  j 


+  Bh  [/]  l  + 


(5.71) 


2  h 


2  h 


=  0 


For  solving  the  considered  system  of  equations  we  need  to  introduce  the  boundary  conditions 
at  each  layer  boundary 

<5'72) 


E,  =E, 

*  p-  *  p+ 


Ez  =Ez 

2  P~  2  P+ 


Hz  =Hz 

2  p-  2  p+ 


(5.73) 

(5.74) 

(5.75) 


From  equation  (5.24)  and  the  following  supstitution  for  coefficients 

kn  1 


Mf  = 


mi  ~kl  p 


r-1  CQSftS 

Mf  =-i — 

m,/ip  ~K 

we  get  the  equation  for  H(i  component  at  the  boundary  between  two  layers 

dE" 

rE  u+z 


h;=m£h:+m: 


dp 


By  discretizing  the  equation  and  by  expressing  the  p  coordinate  as  a  +  ih  we  get: 


h;[n]  =  m^h:[n]+m> 


e:[n]-e:[n-  i] 


(5.76) 


(5.77) 


(5.78) 


(5.79) 


The  boundary  condition  for  the  continuity  of  the  H ^  component  will  give  the  following 
equation 

r  ,  F  i"fiVl-i"[A^-ll  r  n  -F  En\N  +  2]-En\N  +  l] 

M"  Hz  [TV]  +  -  1  J  ^  1 - i  =  M^Hnz  [jV  +  1]  +  Ml  ^ - "  (5-80) 


The  form  that  is  suitable  for  coding  is  equal 
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~  r  -|  |  A/f  |  -  r  n 


Mi 


M 


+  E:[N  +  l]\—L\  +  Ei[N  +  2]\ M  + 


(5.81) 


S:  [  *  ]  } + 6;  [n + 1]  {-m"  }  =  0 


The  boundary  condition  for  Ea  results  in  a  dual  equation: 


El  =  N*En_  +  N: 


and  in  a  dual  discretized-form  equation 


m: 


dp 


(5.82) 


N?  I  ~  r  J  N, 

■+h';[n] 


H 


N: 


H 


h  j  —  h 

&  M{yf  )+e:  [jv+i]{-jyf} = o 

where  abbreviations  represents: 


+h:[n+i]\^\+h:[n+2]< 


Ni 


•+ 


m=- 


kzn  1 


rH  _  :  (°E,EPP 


K  =  j- 


Discretization  of  the  last  two  boundary  conditions  results  in: 

#;[w]=£;[at+i] 

|;[Af]=|;[Af+i] 

Boundary  conditions  at  the  PEC  boundary  are: 

t"[i]  =  0 

t;  id = o 

Boundary  condition  for  the  //-field  together  with  equation  (5.22)  results  in: 


(5.83) 


(5.84) 


(5.85) 


(5.86) 

(5.87) 

(5.88) 

(5.89) 


HnA  1]  =  #[2] 


(5.90) 


For  TMZ  polarization,  the  boundary  conditions  on  the  external  boundary  are  obtained  using 
the  equations  (5. 59. a)  and  (5.59.b): 


Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Metamaterial-Based  Cylinders  Used  for  Invisible  Cloak  Realization 


64 


MHtm[N]  +  M't 


M«sm(a)rC^(kzb)+Mlsm(a)r(y,{kzb)  +  R,Tif\kzb)) 


(5.91) 


NEtt[N]  +  N‘t 


JV?  sin (a)r  (j,(M)+^,H™(i,i))  + JV*  sin(a)y-"C„Hf 1 \kj>) 


(5.92) 


which  results  in 


>  +  i:[N]<  >  +  H:[N]{M^}  +  Rn{-r  sm(a)M*n'n(2\kzb)}+  ^  ^ 

C„  sin(a)M"H<2)(*zZ>)}  =  fn  sm{a)M°Vn{kzb) 

f  ath  1 

+#[*]  [+l;[Jv]{Jv;}+^{-rs>>,(a)Jv/H!,2,(^)}+ 

l  J  l  J  (^*^4) 

C„  Sin(a)jV"Hf>(*2&)}  = sin(a)Np  n(kzb) 


Boundary  conditions  (5.74)  and  (5.75)  results  in: 

I*"  M  =  r  sHa){h(,Kb)+RnK(?{kzbj)  (5.95) 

Hnz  [TV]  =  r  sin(«)C„H^2)(A:zh)  (5.96) 

which  give  us  the  final  two  equations: 

t  sin(a)H[2)(£zZ>)}  =  f"  sin(a  )J„(£z6)  (5.97) 

Hnz  [TV]  +  C„  {-r  sin(«)Hf(A:z6)}  =  0  (5.98) 


For  TEZ  polarization,  the  boundary  conditions  on  the  external  boundary  are  obtained  using  the 
equations  (5. 60. a)  and  (5.60.b): 


M'rH'![N]  +  M‘ 


tz[N]-tz[N-l\ 


m;  sin  (a)fn  (j  „(M)  +  CX2)(M>))  +  m;  sin(a)rRX2\Kb) 


(5.99) 


neJ:[n]+n> 


h:[n]-h:[n-  i] 


TV;  sm{a)j-nRn^\kzb)  +  <  sin(a)y-"  (j^*)  +  Cn^2\kzb)) 


(5.100) 
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which  result  in: 


C,  {-j-  sin(a)M”H«>(^*)j  =  j-  sin(a)M"j.(*,*) 

^[Ar-,]|_^|  +  4-[Af]|^|  +  |;[A,]{^}  +  ^{_rsi„(„)^Hm(t!A)}  +  ^  ^ 
C,  {-j-  sin(a)jv;Hf  >(«)}  =  f  sin (aOJVj'j;^) 


Boundary  conditions  (5.74)  and  (5.75)  result  in: 

t  [#]  =  r  sin(a)RnR^(kzb)  (5.103) 

Hnz  [TV]  =  r  ^(a){in(kzb)+Cn^\kzb))  (5.104) 

which  give  us  the  final  two  equations  in  TEZ  case: 

t  [#]+*,  {-rn  ^m(a)^\kzb)}  =  0  (5.105) 

m  [N]  +  Cn  Sin(a)U[2\kzb)}  =  r  sin(a)J„(*zZ>)  (5.106) 


The  expression  for  the  scattering  width  rroo  and  for  the  total  scattering  width  ctj  are  obtained 
as  follows.  From  the  definition  [17]: 


c t2d  =  lim  2np 

/?— >oo 


E‘ 


(5.107) 


it  follows  (we  substitute  the  expressions  for  the  incident  and  scattered  waves): 


<j2D  =  lim  27 ip 

p — >oo 


oo 

E„  £ 

fl =-00 _ 

E‘e~Ax 


(5.108) 


Using  large  argument  approximation  of  the  Hankel  function  [36] 

Hl2)(*z/?)  =  .pp-/e-^  (5.109) 

]/  nkzp 

we  obtain  the  expression  for  the  scattering  width  cfid- 
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® 2D 


4 

k0  sin(«) 


oo 

z  v* 


For  TEZ  polarization  the  used  definition  for  radar  cross  section  is  [17]: 


cr2D  =  lim  2np 

/?— >00 


(5.110) 


(5.111) 


With  the  large  value  approximation  from  (5.109)  and  with  the  definition  of  incident  and 
scattered  wave  for  TEZ  case  (see  eqs.  (5. 60. a)  and  (5.60.b)),  equation  (5.1 1 1)  is  reduced  to: 


&2D 


4 

k0  sin(«) 


The  total  scattering  width  is  defined  as 


(5.112) 


(5.113) 


where 


represents  the  magnitude  of  the  Poyting  vector  of  the  incident  wave  and  Pl 


represents  the  scattered  power  per  unit  length  (in  the  axial  direction): 


_..|2 


S'  =-L— Lk' 

2?7 


(5.114) 


271 

PL=p\s-f>dt  =  p[ 


2n  2 


nk0p 


iruj) 


(5.115) 


S*=^k'' 
2  rj 


(5.116) 


Using  the  orthogonality  of  complex  exponentials  exp(jn^)  and  after  short  manipulation  we 
get: 


u  n=- oo 


(5.117) 


fn 

n~ — oo 


(5.118) 
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5.2.  RESULTS 

Using  the  developed  program  for  the  oblique  incidence  we  have  tested  the  sensitivity  of 
the  Schurig  cloak  on  the  angle  of  incidence  wave.  The  original  problem  was  a  2D  problem, 
i.e.  it  was  assumed  that  the  incoming  angle  of  the  incident  wave  is  90  degrees.  However,  for 
practical  realizations  it  is  very  important  to  check  the  sensitivity  of  the  considered  cloak  on 
incoming  angle  of  the  incident  wave.  Figure  5.5  illustrates  the  sensitivity  of  the  Schurig 
cloak  (we  have  used  a  10-layer  model  of  the  cloak).  It  can  be  seen  that  the  cloak  has 
scattered- field  reduction  property  up  to  incoming  angle  ~  75°  (or  -105  °).  For  angles  smaller 
than  75°  the  cloaked  cylinder  scatters  more  power  than  the  bare  PEC  cylinder.  That  can  be 
seen  in  Fig.  5.6  where  the  scattered  width  of  the  cloaked  and  bare  PEC  cylinder  is  compared 
for  incident  angle  75°.  The  cross-polarized  scattered  field  is  also  shown  in  Fig.  5.6  (the  cross- 
polarized  scattered  field  is  not  excited  only  for  normal  incidence).  In  Figures  5.7  and  5.8  we 
have  investigated  properties  of  the  TEZ  (Cai)  cloak.  The  parameters  of  the  cloak  are  defined 
by  the  equaton  (3.1.c),  and  we  have  analyzed  a  10-layer  realization  of  the  TEZ  cloak.  It  can  be 
seen  that  the  TMZ  cloak  is  less  sensitive  on  the  angle  of  arival  of  incident  wave  than  the  TEZ 
cloak.  From  Figures  5.7  and  5.8  it  can  be  seen  that  for  TEZ  cloak  for  angles  smaller  than  80° 
the  cloaked  cylinder  scatters  more  power  than  the  bare  PEC  cylinder. 

The  sensitivity  to  angle  of  incidence  can  ba  also  illustrated  by  plotting  the  total  scattered 
width  of  the  PEC  cylinder  with  and  without  the  cloak.  Both  TMZ  and  TEZ  cases  are  considered 
in  Figures  5.9  and  5.10.  It  can  be  seen  that  the  TMZ  (Schurig)  cloak  reduces  the  scattered 
power  for  angles  of  incidence  69°  -  111°  ,  while  the  TEz  (Cai)  cloak  reduces  the  scattered 
power  for  angles  of  incidence  77°  -  103°. 

We  have  also  optimized  the  cloak  dimensions  for  the  angle  of  incidence  that  is  different 
from  the  normal  incidence.  As  an  example  we  have  considered  a  cloak  with  the  reduced 
variation  of  constitutive  parameters  for  TMZ  polarization  (i.e.  the  cloak  which  structure  is 
similar  to  the  Schuring’s  cloak).  In  table  5.1  we  have  compared  the  cloak  parameters 
optimized  for  incidence  angle  90°  and  75°,  respectively.  The  difference  between  structure 
parameters  is  not  too  large,  i.e.  both  structures  have  the  same  level  of  complexity.  The 
comparison  between  scattered  widths  of  these  two  structures  is  given  in  Figure  5.9  (both  of 
them  are  excited  form  the  direction  used  in  the  cloak  design).  It  can  be  seen  that  larger 
reduction  of  scattered  field  can  be  obtianed  for  normal  incidence,  i.e.  when  the  incident  wave 
follows  the  symmetry  of  the  cloak. 
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Figure  5.5.  Normalized  scattering  width  vs.  angle  of  incidence  for  the  Schurig  cloak. 


Figure  5.6.  Comparison  of  the  normalized  scattering  width  of  the  Schurig  cloak  and  of  the 

bare  PEC  cylinder  (incident  angle  is  75°). 
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Figure  5.8.  Comparison  of  the  normalized  scattering  width  of  the  TEZ  cloak  and  of  the  bare 

PEC  cylinder  (incident  angle  is  80°). 
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Figure  5.10.  Total  scattering  width  Gj/2a  vs.  angle  of  incidence  for  the  TEZ  cloak 
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Layer 

No. 

Optimized 
cloak  for 

einc  =  75° 

Optimized 
cloak  for 

Qinc  =  90° 

=3.433; 

=3.841; 

^=1-0 

A>=l-0 

Mr 

Mr 

1 

0.065 

0.021 

2 

0.156 

0.123 

3 

0.255 

0.218 

Table  5.1.  Relevant  constitutive  parameters  of  the  3 -layer  cloaks  optimized 
for  incidence  angles  6inc  =75°  and  6inc  =90°,  respectively. 


Figure  5.11.  Normalized  scattering  width  a2D /2  of  two  cloaks  optimized  for  incidence  angles 
90°  and  75°,  respectively.  Both  cloaks  are  excited  form  the  directions  used  in  the  cloak  design. 
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6  ANALYSIS  OF  UNIAXIAL  MULTILAYER  SPHERICAL 
STRUCTURES  USED  FOR  INVISIBLE  CLOAK 

REALIZATION 
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Analysis  of  uniaxial  multilayer  spherical  structures 

used  for  invisible  cloak  realization 


6.1.  INTRODUCTION 

The  already  presented  cloak  designs  and  optimizations  are  based  on  the  analysis  of  plane 
wave  scattering  by  a  circular  cylinder,  which  is  a  two-dimensional  problem.  Such  approach 
simplifies  the  required  procedure  and  enables  reasonably  accurate  approximation  for  cylinder¬ 
like  objects.  However,  for  scatterers  of  arbitrary  geometry  two-dimensional  approximation  is 
not  general  enough,  which  leads  to  requirement  for  full  three-dimensional  scattering  analysis. 

This  chapter  presents  the  generalization  of  the  analysis  and  optimizations  to  three 
dimensions.  The  analysis  focuses  on  calculating  bistatic  scattering  cross  sections  for  various 
cloak  realizations  and  investigating  the  possibility  of  optimizing  cloak  parameters  to  improve 
the  invisibility  performance.  For  these  purposes,  the  program  for  analyzing  spherical 
structures  was  developed  and  the  optimization  algorithm  was  connected  with  the  extended 
version  of  the  program. 


-\N\ — > 

Incident  plane  wave 


Figure  6.1.  The  cross  section  of  the  analyzed  spherical  structure. 
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6.2.  THEORETICAL  CONSIDERATIONS 

6.2.1.  Vector  eigenvectors 

When  solving  some  differential  equation  such  as  e.g.  wave  equation,  one  of  the  common 
methods  is  to  find  some  appropriate  set  of  linearly  independent  functions  (eigenfunctions)  in 
terms  of  which  it  is  possible  to  express  any  solution  (function).  By  applying  such  method,  to 
obtain  some  particular  solution  (due  to  source  terms  or  boundary  conditions)  one  only  needs 
to  find  the  coefficients  corresponding  to  each  term  of  the  eigenfunction,  i.e.  differential 
equation  could  be  in  principle  be  reduced  to  algebraic  one. 

In  this  section  we  will  show  how  to  construct  the  eigenfunctions  for  vector  wave  equation 
that  is  satisfied  by  the  electric  and  magnetic  field.  Let  C=C(r)  be  some  vector  field  that 
satisfies  the  homogeneous  vector  Helmholtz  equation  (it  represents  e.g.  electric  or  magnetic 
field  or  vector  potential  in  the  absence  of  sources): 

V2C  +  k2C  =  0  (6.1) 

where  k  is  called  wavenumber  (in  general  it  is  a  complex  number). 

As  mentioned  above,  the  idea  is  to  find  some  complete  set  of  vectors  in  terms  of  which  it  is 
possible  to  expand  any  vector  C  that  satisfies  (6.1).  As  a  basis  for  such  expansion  we  first  find 
the  scalar  function  y/  which  satisfies  the  scalar  Helmholtz  equation  (subject  to  some  boundary 
conditions): 


V2y/  +  k2y/  =  0  (6.2) 

Using  scalar  function  <//  and  defining  some  constant  vector  a  it  is  possible  to  derive  three 
sets  of  mutually  linearly  independent  vector  functions  [41]  each  of  which  individually 
satisfies  (6.1.): 


L  =  V !// 

M  =  Vxa^  (6.3) 

N  =  --(VxVxa^)  =  --(VxM) 
k  k 


The  L,  M  and  N  vectors  also  have  the  following  properties  which  arise  due  to  their  definition: 


a)  L  is  irrotational,  i.e.  ( V  x  L)  =  0 

b)  M  and  N  are  solenoidal,  i.e.  V  •  M  =  V  •  N  =  0 

c)  Since  a  is  presumed  constant,  vector  M  could  also  be  written  as:  M  =  L  x  a 

d)  Vector  M  is  perpendicular  to  L,  i.e.  L  •  M  =  0 .  Also,  M  is  perpendicular  to 


]_ 

k 


•(VxN) 


a. 
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If  by  solving  (6.2)  we  find  some  set  of  scalar  functions  {y/n}  which  form  a  complete  set  in  the 
Hilbert  space,  it  can  be  shown  that  any  vector  field  C  that  satisfies  (6.1)  could  be  represented 
by  a  linear  combination  of  L,  M  and  N  vectors  that  correspond  to  each  scalar  t//„  in  the  set,  i.e. 

c=ZKM.+6.N.+c-L.)  «s-4) 


where  a„,  b„  and  c„  are  coefficients  to  be  determined. 

Now  we  turn  our  attention  to  the  electric  and  magnetic  fields  in  region  without  sources  where 
they  are  described  by  the  following  equations  (we  assume  time-harmonic  case  and,  in  first 
case,  the  isotropic  medium): 

VxE  =  -jcojuH 

V  x  H  =  jcoeE  (6.5) 

V-E  =  V-H  =  0 

Since  in  the  source  free  regions  divergence  of  the  electric  and  magnetic  field  is  zero,  it  is 
possible  to  write  them  as  a  curl  of  another  vector  (divergence  of  the  curl  is  always  zero).  If  we 
define  vectors  A  and  F  such  that: 

H  =  —  V  x  A 

**  (6.6) 

E=--VxF 

£ 

for  each  vector  we  can  derive  an  independent  solution  for  both  electric  and  magnetic  field,  by 
inserting  them  into  Maxwell’s  equations  and  applying  Lorenz  condition.  Physically,  vectors 
A  and  F  represent  potentials  that  arise  due  to  electric  and  magnetic  sources,  respectively.  The 
total  field  in  the  region  without  sources  turns  out  to  be  the  superposition  of  the  fields  that 
arise  from  A  and  F : 


E  =  E^  +Ef  = — - — VxVx  A-  —  VxF 

jOJJLIS  £ 

H  =  H,+Hf=-VxA - - — V  x  VxF 


(6.7) 


M 


J  CO jj£ 


Now  we  make  a  particular  choice  for  A  and  F: 


a  y^TM a 
F  =  \j/  TE  •  a 


(6.8) 


where  y/TM  and  y/TE  are  scalar  functions  that  both  satisfy  Helmholtz  equation  (6.2)  and  a  is 
some  constant  vector.  The  expression  (6.7)  could  now  be  rewritten  as: 
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E  =  V x  V x (y/TM  •  a) -- V x  (y/TE  ■  a) 

JCOjUS  s 

H=-Vx(fm’a)-  —  VxVx(^-a) 

jU  J  CO JU8 


(6.9) 


Comparing  the  (6.9)  with  (6.3)  we  immediately  find  that  it  is  possible  to  define  equivalent 
vector  functions  e.g.  Mr£=V x(t//rE  -a)  and  N TM  =  V x  V x ( y/m  ■  a )  in  terms  of  which  it  is 

possible  to  expand  the  electric  and  magnetic  field  in  the  way  as  described  in  (6.4).  Such 
expansion  is  of  the  form: 


E  —  ^(a„Mr£>J1  +  b„NTMn )  —  Er£  +  ErM 

n 

H  =  :  +  bnMTM  )  =  Hr£  +  HrM 

JCOjU  „ 


(6.10) 


We  also  note  that: 

a)  the  expansion  (6.10)  is  a  bit  more  general  than  (6.4)  in  the  sense  that  two  generating 
scalar  functions  (y/TM  and  y/TE)  are  actually  needed  (however  they  only  differ  to  some 
multiplicator). 

b)  in  such  expansion  terms  in  vector  L  are  zero,  which  is  expected  since  divergence  of  L 
by  definition  is  not  generally  zero  unlike  the  divergence  of  electric  and  magnetic  field. 

c)  Since  vector  M  is  perpendicular  to  the  constant  vector  a  we  denote  the  respective  terms 
of  the  field  expansion  in  M  (and  accordingly  their  N  counterparts)  as  transverse  electric 
(TE)  or  transverse  magnetic  (TM)  to  a.  This  means  that  each  field  is  decomposed  to  the  TE 
part  and  TM  part  (generated  by  i//te  and  y/jM ,  respectively). 

To  summarize  the  above  considerations,  we  note  that  the  problem  of  solving  for  electric  and 
magnetic  fields  can  be  formulated  as: 

a)  fixing  the  appropriate  constant  vector  a 

b)  finding  the  complete  sets  of  scalar  functions  y/TE  and  y/TM  by  solving  homogeneous 
scalar  Helmholtz  equation 

c)  expanding  the  field  in  TE  and  TM  parts  (so-called  magnetic  and  electric  multipoles, 
respectively)  and  determining  coefficients  such  as  an  and  bn  by  applying  the  given 
boundary  conditions  or  expanding  the  given  source  functions. 


6.2.2.  Spherical  harmonics 

The  scalar  functions  y/  in  the  expansion  of  the  fields  (6.10)  as  well  as  suitable  choices  of 
vector  a  depend  on  the  choice  of  coordinate  system.  In  this  section  we  describe  procedure  for 
spherical  coordinate  system  ( r,(p,0 )  and  seek  for  expansion  of  the  fields  in  free  space  in  terms 
of  spherical  eigenfunctions,  i.e.  spherical  harmonics.  The  procedure  described  above  requires 
finding  some  fixed  vector  a.  Although  any  vector  a  can  be  chosen  in  principle,  in  spherical 
coordinate  system  it  is  not  possible  to  find  such  vector  that  would  a  priori  generate 
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independent  vectors  M  and  N  across  the  whole  surface  of  the  sphere  -  if  e.g.  we  use  radial 
unit  vector  r ,  we  find  that  it  is  not  constant  across  the  sphere.  Similar  issues  arise  also  in 
more  general  coordinate  systems. 

Nevertheless,  for  spherical  coordinate  systems  in  particular  it  can  be  found  that  from  radial 
vector  it  is  actually  still  possible  to  construct  appropriate  vectors  M  and  N  that  are 
independent  across  the  spherical  surface.  We  try  to  construct  vector  M  tangential  to  sphere 
(which  is  proportional  to  e.g.  electric  or  magnetic  field)  of  the  following  form: 

M  =  V  x.Yu(r)y  (6.11) 

i.e.  we  have  chosen  non-constant  vector  of  the  form  a=u(r)-r ,  where  it(r)  is  some  radial 
function,  yet  undetermined.  After  some  derivation  it  can  be  shown  that  with  selection  u(f)  =  r 
the  function  i//  satisfies  the  scalar  Helmholtz  differential  equation  [41],  [42], 

Now  we  need  to  find  the  complete  set  of  functions  i//  that  satisfy  scalar  Helmholtz  equation 
(6.14c)  in  spherical  coordinates.  The  typical  procedure  is  the  method  of  separation  of 
variables.  We  assume  that  the  function  i//=  i//(r)  can  be  written  as  a  product  of  three  functions 
of  r-,  9-  and  cp-  coordinates,  respectively: 

Y  =  f(r)-g(6)-h(cp)  (6.12) 


This  form  of »//  is  then  inserted  in  Helmholtz  equation  in  spherical  coordinates,  and  after  some 
manipulations  one  obtains  three  independent  equations  for  functions/  g  and  h : 

a)  ID  Helmholtz  equation 


d2h  2 


d(p 2 


+  m  h  =  0 


b)  Legendre  differential  equation 
1  d 


+ 

n{n  + 1)  - 

m  x 

2 

89  _ 

v  sin  9  j 

sin  6  dO  \ 

c)  Spherical  Bessel  differential  equation 


8  =  0 


8  f  2  df 

- 2+ - ~  + 

dr  r  dr 


2  n(n  + 1) 


k- 


f  =  o 


(6.13a) 


(6.13b) 


(6.13c) 


Here  m 2  and  n(n+ 1)  are  separation  constants,  while  m,n  are  integers.  The  solutions  of  (6.13a) 
are  harmonic  functions  h(<p)  =  e~jm<p ,  while  the  solutions  (6.13b)  are  Legendre  functions  of  the 

first  and  second  kind  g{9)  =  P™ (cos#)  and  g(9)  =  Q™ (cos 9)  (the  latter  ones  give  the  non¬ 
physical  solution  due  to  singularity  at  9  =  0°  and  9  =  180°).  Finally,  the  solutions  of  (6.13c) 
are  spherical  Bessel  and  Hankel  functions  of  the  first  and  second  kind  / (r)  =  ( kr ) 
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and / (r)  =  h(2)  (kr) .  We  choose  solutions  in  accordance  to  physical  reality  of  the  problem  in 

question  (outgoing  wave  in  free  space,  no  singularities  on  the  spherical  surface)  and  write  the 
normalized  solution  for  y/  as: 

Vm,  =  /M  •  g(0)  ■  m  =  €  (V)  ■  p;  (cos  0)  ■  e-P*  (6. 14) 

Note  that  actual  solution  is  actually  multiple  of  y/mn.  If  one  needs  to  take  into  account  standing 
waves  in  some  spherical  layer,  then  spherical  Bessel  functions  j„(kr)  should  be  also 
considered. 

Spherical  vector  eigenfunctions  L,  M  and  N  can  be  now  constructed  from  each  y/mn,  i.e.  by 
inserting  (6.14)  into  (6.11)  and  (6.3).  Therefore,  we  have  by  now  obtained  proper 
mathematical  tool  to  expand  any  vector  function  in  spherical  coordinates  as  a  sum  of  vector 
eigenfunctions.  The  remaining  task  is  to  apply  the  developed  tool  to  potential  and  field 
equations  and  to  obtain  required  field  expressions  in  terms  of  vector  eigenfunctions. 

In  accordance  with  discussions  described  above  we  expand  potentials  A  and  F  as: 

„  (6.15) 

Y  =  yrTF.y  =  ^AnxVmn-rx 

m,n 


where  Amn  and  Bmn  are  coefficients  to  be  determined. 

We  note  that  the  potentials  contain  only  the  radial  component.  The  final  expansion  of  the 
electric  and  magnetic  field  is  then  obtained  by  inserting  (6.15)  into  (6.9)  which  is  the 
expression  for  electric  and  magnetic  field  in  terms  radial  components  of  potentials  that 
provide  decomposition  to  TE  and  TM  (to  r-coordinate)  parts.  The  expression  is  of  the  form 
(6.10)  only  with  double  index  mn.  The  coefficients  in  the  expansion  are  then  obtained  by 
applying  boundary  conditions  or  by  expanding  the  (known)  source  terms  in  spherical 
harmonics  and  finding  the  field  coefficients  that  correspond  to  each  source  term  (i.e. 
differential  equation  is  then  reduced  to  algebraic  one). 

As  a  final  remark,  the  radial  dependence  in  (6.14)  is  often  incorporated  in  modified  form  of 
spherical  Hankel  functions  (so-called  Schelkunoff  form): 


H(2)  ( kr )  =  kr  ■  h(2)  ( kr  ) 


(6.16) 


The  connection  between  spherical  Hankel  functions  in  Schelkunoff  and  ordinary  form  H{2] 
and  h(2)  with  the  regular  (cylindrical)  functions  H(2)  is  given  by: 


H(2)  (At*)  =  kr  •  h{2)  (kr)  =  kr  ■ 


(6.17) 
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The  same  relation  as  (6.17)  holds  also  for  other  types  of  Bessel/Hankel  functions.  The 
Bessel/Hankel  functions  of  the  Schelkunoff  form  are  solutions  to  the  slightly  modified 
equation  compared  to  (6.13c): 


d2f 
dr 2 


+ 


k2  n(n  + 1) 


/  =  o 


(6.18) 


Note  that  the  expansions  of  the  type  (6.15)  are  only  slightly  changed  by  employing  the 
Schelkunoff  form  of  Hankel  functions  instead  of  the  ordinary  ones  (only  the  coefficients  such 
as  Amn  or  Bmn  are  different  and  r  is  incorporated  in  Hankel  functions). 


6.2.3.  Modifications  for  anisotropic  structures 

In  this  section  we  consider  the  extension  of  the  above  described  method  of  field  expansion 
to  the  anisotropic  structures.  In  particular,  we  consider  uniaxial  structures  where  electric 
permittivity  and  magnetic  permeability  are  given  by  the  tensors  of  the  form: 


£  =  £n 


£r 

0 


0  0 

£,  0 


0  0  e. 


A  =  A0 


jur  0  0 

0  *  0 
. 0  0  A. 


(6.19a) 


(6.19b) 


which  means  that  the  radial  (r-)  and  transverse  (<p-  and  6-)  components  of  the  fields  “see” 
different  properties  of  the  medium  when  propagating  through  such  media.  The  curl  Maxwell 
equations  in  the  region  without  sources  are  now  of  the  form  similar  to  (6.5): 


VxE  =  -jojjLiW 
V  x  H  =  ja>£ E 


(6.20) 


The  idea  is  to  expand  electric  in  magnetic  fields  in  spherical  eigenvectors  according  to 
procedure  described  in  previous  sections.  To  make  such  expansion,  one  needs  to  obtain 
proper  differential  equation  for  the  anisotropic  case.  To  make  dealing  with  tensors  more 

convenient,  we  make  use  of  the  flux  density  vectors  B  =  //H  and  D  =  sE  rather  than  fields  E 
and  H,  and  thus  rewrite  the  curl  equations  as: 


Vx 

Vx 


s  D 


//  B 


=  -j<*> B 


=  jcoD 


(6.21) 
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In  accordance  to  the  procedure  previously  described,  we  now  choose  potentials  with  radial 
components  only,  by  which  we  construct  independent  expressions  for  the  TE  (from  Fr)  and 
TM  (from  Ar)  fields.  The  potentials  (6.15)  are  therefore  written  as: 


A  =  Ytm  ■rr  =  Ar  -r 
F  =  i//rE  -rr  =  F-r 


(6.22) 


while  with  such  choice  of  potentials  the  TE  and  TM  flux  density  vectors  are  then  given 
corresponding  to  (6.6)  as 


Btm  =  Vx  A  =  V x(Ar  -r) 

D tr  =-VxF  =  -Vx(Fr  -r) 


(6.23) 


Now  we  substitute  (6.23)  into  (6.21)  and  obtain  respective  TM  and  TE  counterparts  of  (6.23): 


B  7  =  — V  x 


D 


TM 


JO) 

J_ 

jO) 


C=- 1 

£  -(Vxf^f) 


V 


(=- 1 


Vx 


A  •  ( V  x  Arr ) 


(6.24) 


Using  (6.21)  we  can  choose  to  establish  the  following  relations  between  TM  and  TE  parts  of 
D  and  B  and  thus  obtain  decoupled  expressions  in  terms  of  A  and  F,  respectively: 


Vx 

Vx 


£  D 

ji  B 


TM 


TE 


—  jcoBm 
=  jo)J)TE 


(6.25) 


The  equations  (6.25)  are  now  expanded  using  (6.24)  and  (6.23)  and,  since  radial  component 
of  the  left  sides  is  zero  (due  to  definition  of  TM  and  TE  cases),  after  some  algebraic 
manipulations  we  arrive  to  the  differential  equations  for  Ar  and  Fr  (i.e.  for  ijtm  and  i//te)'- 


£r  d2Ar  1  d 

— — - -H - 

st  dr2  r 2  sin  6  86 

Mr  ,  1  d 

/ut  dr2  r 2  sin  6  86 


sin^ 


sin^ 


a i 

86 

8F]_ 

86 


1  82Ar 

H - 

r2  sin2  6  8cp 2 
1  82Fr 

H - 

r2  sin2  6  8(p 2 


+  k2  fj.tsrAr  =  ° 
+  k2jir£tFr=0 


(6.26a) 

(6.26b) 


where  k2  =  or  jins() . 

By  comparison  with  Helmholtz  equation  in  spherical  coordinates  (6.13)  we  note  that  the 
expressions  (6.26)  can  be  easily  shown  (by  doing  some  rewriting)  to  be  actually  equivalent, 
except  for  the  factors  sj st  and  ji,J ji,  (so-called  anisotropy  ratios)  in  the  term  of  the  second 
radial  derivative.  We  also  note  that  TE  and  TM  functions  do  not  satisfy  the  same  equation  in 
terms  of  medium  parameters,  as  was  the  case  in  isotropic  media. 
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Using  the  separation  of  variables  method,  one  obtains  the  same  9-  and  <p-  dependence  as  in 
(6.13a)  and  (6.13b),  while  for  r-dependence  the  following  equations  for  radial  dependences  of 
Ar  and  F,  are  obtained: 


d2Ar 
dr 2 


st  n(n  + 1) 


(6.27a) 


d2Fr 
dr 2 


+ 


ko  stHt 


Mj  njn  + 1) 
Mr  r 2 


(6.27b) 


These  equations  turn  to  be  spherical  Bessel  differential  equations  of  Schelkunoff  form  (6.18). 
The  solutions  are  then  spherical  Bessel/Hankel  functions.  For  both  cases  the  argument  is: 

K^Mftr  =  ktr  (6.28) 


On  the  other  hand,  the  orders  of  Bessel/Hankel  function  is  different  for  each  case: 


a)  for  TM  case: 


vi  = 


n(n  +  l)  —  +  — 
s  4 


2 


b)  for  TE  case: 


n(n  + 1)^  +  -  2 
Mr  4_ 


\_ 

2 


(6.29a) 


(6.29b) 


We  also  note  that  the  resulting  orders  are  not  integers.  Now  we  can  construct  the  set  of 
eigenfunctions  in  terms  of  which  expand  radial  vector  potentials  are  expanded  according  to 
procedure  described  for  isotropic  case.  The  expansions  are  of  the  form: 

A  =  Vtm?  =  {k,ryP;(coS0)-e-" 

(6.30) 

K  =  Vte r  =  £  (k.’-yC  (cosO)  ■  e** 

m,n 

The  Hankel  functions  of  the  second  kind  represent  waves  traveling  outside  the  sphere.  For 
fields  expansions  inside  the  sphere  one  only  needs  to  add  terms  with  Bessel  functions  of  the 
first  kind,  Jv  and  ,Jv  ,  which  represent  the  standing  waves. 

The  remaining  task  is  to  construct  vector  eigenfunctions  M  and  N  and  expand  the  electric 
and  magnetic  field  in  terms  of  them.  Since  we  are  dealing  with  anisotropic  media  the 
intermediate  step  consists  of  expressing  the  flux  density  vectors  B  and  D  according  to  (6.23) 
and  (6.24)  by  which  the  actual  fields  follow  due  to  properties  of  TE  and  TM  cases  as: 
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E  = — D 

i^Te  ** TE 


it  _  _  n 

nrM  D7¥ 


(6.31) 


The  rest  of  the  procedure  is  the  same  as  in  the  isotropic  case  described  in  the  previous  section, 
while  the  final  expansion  is  of  the  form  as  given  in  (6.10),  only  with  double  indexes  ( mn ). 


6.2.  RESULTS  OF  THE  CLOAK  ANALYSIS 

6.2.1.  Pendry  cloak  design 

For  analysis  of  cloaking  of  three-dimensional  objects,  cloaking  of  a  PEC  sphere  of  0.75k 
radius  (i.e.  a=2.645  cm  at  the  frequency  of  8.5  GHz)  was  considered  and  few  configurations 
were  compared.  The  achieved  invisibility  is  expressed  in  terms  of  bistatic  scattering  cross 
section  030  (normalized  to  the  cross-section  of  the  PEC  sphere)  [17],  In  the  first  configuration, 
the  cloak  constitutive  parameters  have  been  calculated  by  using  coordinate  transformations 
given  in  [9]  (so-called  Pendry  design)  which  required  radial  variations  of  radial  components 
of  the  permittivity  and  permeability  tensor,  respectively,  while  keeping  transverse 
components  constant: 


G  =Mr 

£n  =  /C 

£f=^ 


b  ( r—a^ 


b-a 

b 

b-a 

b 

b-a 


V  r  ) 


(6.32) 


These  variations  were  approximated  by  piecewise  constant  functions  representing  the  number 
of  cloak  layers.  In  Fig.  6.2  comparison  of  the  obtained  bistatic  cross  section  for  the  uncloaked 
PEC  sphere  and  the  one  cloaked  with  Pendry  design  is  shown.  The  influence  of  the  cloak 
thickness  was  also  taken  into  account.  For  thick  cloak  design  (outer  radius  of  2a)  the  forward 
scattered  field  (0°  direction)  is  even  larger  than  the  one  of  the  PEC  sphere  itself  which 
actually  makes  it  unsuitable  for  cloaking  purposes.  On  the  other  hand,  for  thin  cloak  design 
(outer  radius  of  1.06a)  both  forward  (shadow)  and  backward  (reflection)  scattered  field  are 
significantly  reduced,  since  thin  cloaks  occupy  smaller  region  of  space.  Therefore  thin  cloaks 
are  more  interesting  for  further  study. 
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Figure  6.2  Normalized  bistatic  scattering  cross  section  of  PEC  sphere  and  two  versions  of 
Pendry  cloak  design 


6.2.2.  Engheta  cloak 

Another  approach  to  cloaking  is  based  on  expanding  the  scattered  field  in  electric  and 
magnetic  multipoles  and  reducing  the  magnitude  of  several  relevant  multipoles  [43]  via 
optimizing  the  cloak  constitutive  parameters  and  required  thickness  (so-called  Engheta 
cloak).  The  comparison  between  performance  of  the  Engheta  cloak  proposed  in  [43]  and  the 
Pendry-designed  thin  cloak  (5  layers)  is  given  in  Fig.  6.3.  For  both  cloaks  the  outer  radius  is 
1 ,06a.  It  should  be  noted  that  although  Pendry  cloak  design  provides  much  better  invisibility 
performance  (Fig.  6.3),  Engheta  cloak  is  one-layered  and  isotropic  [43],  which  results  in 
much  simpler  design  in  practice. 


Figure  6.3  Normalized  bistatic  scattering  cross  section  of  PEC  sphere,  Engheta  cloak  and  the 
Pendry  (thin)  design 
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6.2.3.  Optimized  cloak 

As  suggested  in  [43],  the  technique  described  above  could  be  improved  by  introducing 
more  cloak  layers  and  allowing  material  anisotropy.  Such  extension  was  carried  out  in  this 
paper.  By  using  the  developed  PSO  algorithm  we  have  considered  cases  with  various 
numbers  of  anisotropic  layers  and  have  optimized  the  permittivity  and  permeability  tensor 
components  in  order  to  achieve  minimum  scattering  cross  section.  The  overall  cloak  thickness 
was  kept  consistent  to  previous  configurations  (i.e.  outer  radius  of  1 .06 a).  It  was  found  that 
with  only  two  layers  it  is  possible  to  find  such  combination  of  parameters  (Table  6.1)  that 
enables  significant  scattering  cross  section  reduction  for  both  forward  and  backward 
direction.  In  comparison  with  Pendry  cloak  design  (5  layers)  given  in  Fig.  6.4,  forward 
scattering  has  been  even  more  reduced  while  the  practical  realization  of  2-layers  optimized 
cloak  should  be  simpler  (although  anisotropic  dielectric  and  magnetic  materials  are  required). 
For  illustration,  in  Fig.  6.5  the  field  plots  around  the  structures  for  the  two  best  obtained  cloak 
configurations  as  well  as  PEC  sphere  are  given  in  terms  of  total  field  (i.e.  incident+scattered 
field),  and  scattered  field  itself. 


Layer  No. 

Outer  radius  [cm] 

Relative 

permittivity 

Relative 

permeability 

Transverse 

Radial 

Transverse 

Radial 

1 

2.7328 

3.9223 

3.0994 

2.3165 

1.8378 

2 

2.8006 

3.9223 

0.9086 

2.3165 

0.8540 

Table  6.1  Constitutive  parameters  of  the  optimized  cloak 


Figure  6.4  Normalized  bistatic  scattering  cross  sections  of  PEC  sphere,  Pendry  (thin)  design 
and  optimized  cloak. 
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A2.  PEC  sphere 


B2.  Pendry  cloak  design  (thin) 


Cl.  Optimized  cloak 


C2.  Optimized  cloak 


Figure  6.5.  Plots  of  total  and  scattered  field  distribution  in  the  vicinity  of  a  PEC  sphere  (A)  and  two 
cloak  configurations  -  five  layers  Pendry  design  (B)  and  two-layers  optimized  design  (C).  The  plane 
wave  is  impinging  from  top.  The  scale  is  normalized  for  all  cases. 
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Conclusions 


During  the  12  months  assigned  to  this  project  the  software  for  full-wave  analysis  of 
uniaxial  multilayer  cylinders  used  for  invisible  cloak  realization  has  been  developed.  The 
analysis  approach  is  based  on  the  modal  approach,  i.e.  distribution  of  the  electromagnetic 
field  in  radial  direction  is  described  with  a  modified  Bessel  differential  equation  that  takes 
into  account  the  radial  anisotropy  of  the  proposed  structures.  Such  an  approach  is  suitable  for 
2D  problems,  i.e.  when  the  incident  wave  is  propagating  in  the  perpendicular  direction  to  the 
structure  axis.  For  the  oblique  incidence  case  we  have  developed  more  general  analysis 
approach  in  which  distribution  of  the  electromagnetic  field  in  radial  direction  is  determined 
using  finite-difference  discretization  of  the  appropriate  partial  differential  equations.  In  order 
to  characterize  the  level  of  “invisibility”  achieved  by  metamaterials,  we  have  calculated  the 
total  scattering  width,  as  well  as  the  angular  variations  of  scattering  width  in  the  azimuthal 
plane.  The  scattering  properties  of  a  PEC  cylinder  are  also  given  as  a  referent  case  necessary 
to  calculate  the  total  scattering  width  reduction. 

In  this  report  we  have  shown  that  the  achieved  total  scattering  width  reduction  is  only 
around  three  for  structures  proposed  by  Pendry  and  Schurig  with  anisotropy  in  radial 
direction  only.  The  reason  for  such  a  small  cloak  gain  lies  in  the  reflection  of  the  incident 
wave  from  the  cloak  surface  due  to  the  impedance  mismatch.  Furthermore,  it  is  shown  that 
the  needed  continuous  radial  variation  of  constitutive  parameters  can  be  successfully 
approximated  with  about  6  layers  of  constant  permittivity  and  permeability,  meaning  there  is 
hardly  any  improvement  of  the  achieved  invisibility  for  structures  with  more  layers.  The 
needed  number  of  layers  with  constant  constitutive  parameters  however  strongly  depends  on 
the  radius  of  the  object  that  we  would  like  to  cloak.  Furthermore,  it  was  shown  that  the  cloaks 
that  are  designed  for  normal  incidence  can  successfully  work  for  incident  angles  in  the  range 
75°  -  105°.  If  the  incident  angle  is  different  from  90°,  our  analysis  approach  enables  us  to 
make  a  cloak  design  for  a  particular  angle  of  incidence  simply  by  connecting  the  developed 
program  with  the  optimization  routine. 

The  dispersion  of  split-ring-resonator-based  structures  and  thin-wire -based  structures  were 
taken  into  account  by  Lorentz  and  Drude  models,  respectively.  It  was  found  that  the  achieved 
invisibility  is,  due  to  the  dispersion,  a  narrow-band  phenomenon.  For  the  TMZ  cloak  (i.e.  for 
the  cloak  that  was  experimentally  made  by  Schurig  et  al .)  the  invisibility  bandwidth  (relative 
frequency  range  where  the  total  scattering  width  reduction  is  larger  than  1)  is  about  0.24%, 
while  for  the  TEZ  cloak  (i.e.  the  cloak  that  can  be  made  using  periodic  wire  media)  the 
invisibility  bandwidth  is  a  bit  larger  and  its  value  is  about  2.4%.  This  is  explained  by  the  fact 
that  the  thin-wire-based  structures  are  less  sensitive  to  frequency  changes  than  the  split-ring- 
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resonator-based  structures.  All  the  results  of  the  presented  analysis  show  that  metamaterial- 
based  cloaks  are  still  in  their  infant  phase,  while  practical  engineering  realizations  are  yet  to 
be  seen. 

We  have  also  merged  the  developed  analysis  program  with  the  global  optimization  routine 
based  on  the  PSO  algorithm.  By  this  we  have  developed  a  software  platform  that  we  have 
used  to  provide  answers  to  two  questions:  is  it  possible  to  realize  a  cylindrical  cloak  that  is 
entirely  invisible  (at  least  at  the  central  frequency)  only  with  the  radial  component  of 
permeability  tensor  being  coordinate  dependent,  and  is  it  possible  to  enlarge  the  bandwidth  of 
such  a  cloak?  The  obtained  results  show  that  it  is  indeed  possible  to  realize  the  “perfect” 
cylindrical  cloak,  but  at  the  expense  of  bandwidth  that  is  reduced  to  just  about  0.08%.  The 
discussion  on  the  influence  of  the  number  of  cloak  layers  shows  that  with  only  three  layers  it 
is  possible  to  obtain  a  better  cloak  performance  than  the  Schurig  cloak  (containing  ten  layers), 
and  that  with  ten  or  more  layers  one  can  design  a  cloak  that  is  practically  invisible.  However, 
these  results  strongly  depends  on  the  cloaked  object  radius,  and  more  layers  are  needed  for 
structures  with  larger  radius  (on  the  other  side,  for  structures  with  small  radius  two-layer 
cloaks  give  very  good  results).  Furthermore,  it  is  shown  that  the  bandwidth  of  the  Schurig 
cloak  can  indeed  be  increased  by  up  to  a  factor  of  2.5,  by  optimizing  the  cloak  layer 
parameters,  while  maintaining  the  same  minimal  value  of  the  total  scattered  width.  Finally,  it 
is  demonstrated  that  losses  in  the  cloak  increase  the  total  scattering  width,  or  in  other  words, 
reduce  the  invisibility  gain.  However,  for  values  of  around  y=  1 0'3/0  the  invisibility  gain  of 

the  optimized  cloak  is  still  roughly  equal  to  the  invisibility  of  the  Schurig  cloak  with  no 
losses. 

Finally,  we  have  extended  the  developed  program  to  uniaxial  spherical  cloaks.  The  three- 
dimensional  analysis  is  in  principle  more  suitable  for  arbitrary  geometries  of  the  object.  We 
have  compared  the  properties  of  cylindrical  and  spherical  cloaks,  and  the  largest  difference  is 
in  the  thickness  of  the  optimum  cloak.  While  in  the  cylindrical  case  the  typical  ratio  between 
the  outer  and  inner  radius  of  the  cloak  is  around  2,  in  the  spherical  case  this  ratio  is  much 
smaller,  typically  around  1.1.  By  merging  the  developed  program  for  analyzing  multilayered 
spherical  structures  with  the  PSO  global  optimization  routine  it  was  found  that  it  is  possible  to 
realize  a  cloak  with  only  two  thin  anisotropic  layers,  which  significantly  reduces  scattered 
field  in  comparison  to  other  cloak  realizations  given  in  literature. 

As  a  short  summary,  the  realized  outcomes  of  the  project  are: 

•  We  have  developed  the  program  “UniaxCloak”  that  analyzes  cylindrical  and  spherical 
cloaks  made  from  the  uniaxial  materials.  The  program  can  fully  characterize  uniaxial 
cloaks,  i.e.  it  calculates  scattered  width  and  total  scattered  width  in  a  frequency  range 
of  interest.  In  the  cylindrical  case,  the  incident  wave  can  have  arbitrary  angle  of 
incidence. 

•  We  have  merged  the  program  “UniaxCloak”  with  the  global  optimization  program 
based  on  the  particle  swarm  optimization  (PSO)  method.  By  this  we  can  determine  the 
optimum  cloak  design  since  the  analysis  methods  based  on  the  transformation 
electromagnetics  do  not  take  all  the  effects  into  account. 

•  We  have  made  a  detailed  investigation  of  properties  of  cloaks  that  are  made  from 
uniaxial  materials.  In  practice,  such  cloaks  can  be  made  from  metamaterials,  for 
example  from  split-ring  periodic  structures  or  periodic  wire  structures. 
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